This supplementary R Markdown document provides the full analytical
workflow and visual outputs referenced in the main manuscript “The
legacy of past heatwaves on ‘off-host’ parasite stages: Reduced
infection risk and costs in the Daphnia-Pasteuria system”. All analyses
were performed in R 4.4.2 (2024-10-31). The dataset and scripts are
provided to ensure full reproducibility.
This section outlines the steps taken to clean, transform and structure our dataset for subsequent analysis. Our dataset comprises various measurements of a sample of biological entities, including spore counts, infection statuses and life cycle metrics such as clutch size, offspring numbers and developmental stages. However, the raw dataset contained missing values that needed to be addressed prior to analysis.
Below is a detailed explanation of each column in the dataset. This breakdown is intended to provide clarity on the variable names, units, and biological significance.
ID: A unique identifier for each sample, used to distinguish individual observations.
Temperature: The experimental temperature (in °C) at which each sample was maintained. The possible values are 20, 30 and 40 °C, which were used to examine the effects of temperature on biological parameters.
Isolate: Identifier for the specific parasite genotype used in the experiment (e.g. C1, C14, C18, C24). Each genotype represents a distinct biological strain or lineage of the parasite.
Clutch: The number of clutch events recorded for each sample, representing reproductive events or cycles.
Nb_offspring: The total number of offspring produced by each sample across clutch events.
FirstClutchAge: The age (in days) at which the first clutch was recorded.
LastClutchAge: The age (in days) at which the last clutch was recorded.
DeathAge: Age (in days) at death, providing a measure of lifespan for each sample.
Accidental_death: Indicates whether the death was accidental. Missing values suggest non-accidental or unspecified causes.
Volume_counted: The volume of the sample (in microlitres) counted under observation and used to standardise spore counts.
InfectedStatus: The infection status of each sample, with possible values including ‘infected’, ‘not_inf’ (not infected) or ‘lost’ for missing data.
Number of mature spores: The raw count of mature spores in the sample.
Number of immature spores (imm):
The raw count of immature spores representing early-stage spore
development.
Nb of white spores (whites): The
number of white spores, which are mature but may appear white instead of
red at the microscope.
cauliflowers: A categorical indicator (yes/no) of whether cauliflower-like structures were observed in the sample.
grapes: A categorical indicator (yes/no) of the presence of grape-like structures, which are another morphological feature observed in spore samples.
mature_load: The calculated load of mature spores, standardised based on ‘Volume_counted’. This represents the density of mature spores in the sample volume.
imm_load: The load of immature spores, calculated similarly to ‘mature_load’ for standardised comparison.
ratio: The ratio of mature spores to the total spore load (mature and immature). This metric provides an indication of spore maturation and infection severity.
Batch: A grouping variable created based on ID, representing experimental batches. This factor enables batch effects to be controlled in the analysis.
group: A combination of isolation and temperature conditions (e.g. C1840) used to categorise samples in specific experimental settings.
Variable Conversion: Several columns required type conversion to facilitate analysis. For example, ‘Volume_counted’ and ‘Nb of white spores’ were converted to numeric values. During this process, some ‘NA’ values were introduced due to incompatible data entries.
Factorization: The categorical variables ‘cauliflowers’, ‘grapes’ and ‘group’ were converted into factor variables. Additionally, the ‘Isolate’ and ‘Temperature’ variables were explicitly factored to ensure clear categorical representation.
Data Subsetting and Exclusion: Observations with specific conditions, such as an ID of 50, were excluded to maintain data integrity.
Calculation of Derived Variables:
Grouping and Batching: Observations were assigned to batches for experimental grouping based on their ID. This enables batch effects to be tracked and controlled in downstream analysis.
## ID Temperature Isolate Clutch Nb_offspring
## Min. : 1.0 20:129 C14: 92 Min. : 1.000 Min. : 0.00
## 1st Qu.: 169.0 30:125 C1 : 87 1st Qu.: 1.000 1st Qu.: 0.00
## Median : 627.0 40:143 C18:116 Median : 5.000 Median : 10.00
## Mean : 517.2 C24:102 Mean : 5.666 Mean : 27.48
## 3rd Qu.: 903.0 3rd Qu.:10.000 3rd Qu.: 57.25
## Max. :1048.0 Max. :16.000 Max. :111.00
## NA's :101 NA's :1
## FirstClutchAge LastClutchAge DeathAge Accidental_death
## Min. : 1.00 Min. :10.00 Min. :15.00 Min. :1
## 1st Qu.:14.00 1st Qu.:21.00 1st Qu.:39.00 1st Qu.:1
## Median :18.00 Median :42.00 Median :52.00 Median :1
## Mean :20.41 Mean :42.49 Mean :49.26 Mean :1
## 3rd Qu.:24.00 3rd Qu.:67.00 3rd Qu.:67.00 3rd Qu.:1
## Max. :67.00 Max. :67.00 Max. :67.00 Max. :1
## NA's :102 NA's :107 NA's :391
## Volume_counted InfectedStatus Number of mature spores
## Min. : 150.0 Length:397 Length:397
## 1st Qu.: 150.0 Class :character Class :character
## Median : 150.0 Mode :character Mode :character
## Mean : 274.9
## 3rd Qu.: 450.0
## Max. :1200.0
## NA's :12
## Number of imature spores Nb of white spores cauliflowers
## Length:397 Length:397 counted-after-freezer: 7
## Class :character Class :character no :270
## Mode :character Mode :character yes : 98
## NA's : 22
##
##
##
## grapes whites group statusDeath
## counted-after-freezer: 7 Min. : 0.000 C1840 : 43 Min. :0.00
## no :266 1st Qu.: 0.000 C2420 : 39 1st Qu.:0.00
## yes :101 Median : 0.000 C1440 : 38 Median :1.00
## NA's : 23 Mean : 6.753 C1830 : 38 Mean :0.67
## 3rd Qu.: 1.000 C1820 : 35 3rd Qu.:1.00
## Max. :348.000 C2440 : 34 Max. :1.00
## (Other):170
## Number of immature spores Never_repro Batch imm
## Length:397 Min. :0.0000 11 : 12 Min. : 0.00
## Class :character 1st Qu.:0.0000 6 : 11 1st Qu.: 0.00
## Mode :character Median :0.0000 28 : 11 Median : 0.00
## Mean :0.2544 4 : 10 Mean : 16.58
## 3rd Qu.:1.0000 8 : 10 3rd Qu.: 18.00
## Max. :1.0000 14 : 10 Max. :233.00
## (Other):333
## mat mature_load imm_load ratio
## Min. : 0.00 Min. : 0 Min. : 0 Min. :0.0000
## 1st Qu.: 0.00 1st Qu.: 0 1st Qu.: 0 1st Qu.:0.7821
## Median : 0.00 Median : 0 Median : 0 Median :0.8734
## Mean : 96.77 Mean : 2357630 Mean : 382870 Mean :0.8179
## 3rd Qu.:198.00 3rd Qu.: 4207500 3rd Qu.: 360000 3rd Qu.:0.9414
## Max. :899.00 Max. :20227500 Max. :4230000 Max. :1.0000
## NA's :12 NA's :12 NA's :250
## repro_rate stade_max
## Min. :0.0000 None :235
## 1st Qu.:0.0000 Cauliflower : 14
## Median :0.2440 Grapes : 1
## Mean :0.4607 Immature : 2
## 3rd Qu.:0.9254 Mature : 68
## Max. :1.6567 2nd_Cauliflower: 77
## NA's :1
Objective: primary goal is to determine the
effect of genotype (Isolate) and temperature conditions
(Temperature) on each of the measured traits. The aim is to
assess the individual and interactive effects of these variables while
accounting for the random effects associated with batch
variability.
Model Selection:
Results and Interpretation: The final model outputs are summarised, enabling us to interpret the fixed effects of temperature on infection likelihood while controlling for batch variability.
This section examines the impact of temperature and the parasite’s genotype on infection rates. Specifically, we use generalised linear mixed models to examine how these factors, along with random batch effects, contribute to the probability of infection.
crushed$InfectedStatus = factor(crushed$InfectedStatus, levels=c("not_inf","infected"))
# Step 1: Evaluating the importance of random effect (Batch)
Inf5 = glmmTMB(data = crushed, InfectedStatus ~ Isolate * Temperature, family = "binomial", REML = TRUE)
Inf5R = glmmTMB(data = crushed, InfectedStatus ~ Isolate * Temperature + (1 | Batch), family = "binomial", REML = TRUE)
# Comparing models with and without random effect using AICc
AICc(Inf5, Inf5R)
## df AICc
## Inf5 12 524.4920
## Inf5R 13 520.4055
# Step 2: Fixed effects selection
Inf0 = glmmTMB(data = crushed, InfectedStatus ~ 1 + (1 | Batch), family = "binomial", REML = FALSE)
Inf1 = glmmTMB(data = crushed, InfectedStatus ~ Isolate + (1 | Batch), family = "binomial", REML = FALSE)
Inf2 = glmmTMB(data = crushed, InfectedStatus ~ Temperature + (1 | Batch), family = "binomial", REML = FALSE)
Inf3 = glmmTMB(data = crushed, InfectedStatus ~ Isolate + Temperature + (1 | Batch), family = "binomial", REML = FALSE)
Inf4 = glmmTMB(data = crushed, InfectedStatus ~ Isolate * Temperature + (1 | Batch), family = "binomial", REML = FALSE)
# AICc comparison for fixed effects models
AICc(Inf0, Inf1, Inf2, Inf3, Inf4)
## df AICc
## Inf0 2 528.0016
## Inf1 5 531.1036
## Inf2 4 510.8766
## Inf3 7 514.1415
## Inf4 13 515.6191
# Model summaries for selected models
tab_model(Inf2, Inf3, Inf4, Inf1)
| InfectedStatus | InfectedStatus | InfectedStatus | InfectedStatus | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| (Intercept) | 0.76 | 0.58 – 0.98 | 0.037 | 0.63 | 0.39 – 1.02 | 0.062 | 0.63 | 0.39 – 1.03 | 0.067 | 0.60 | 0.38 – 0.95 | 0.029 |
| Temperature [linear] | 0.42 | 0.29 – 0.62 | <0.001 | 0.42 | 0.29 – 0.62 | <0.001 | 0.42 | 0.19 – 0.95 | 0.037 | |||
| Temperature [quadratic] | 1.03 | 0.71 – 1.50 | 0.866 | 1.05 | 0.72 – 1.53 | 0.801 | 0.63 | 0.29 – 1.40 | 0.262 | |||
| Isolate [C1] | 1.65 | 0.86 – 3.16 | 0.132 | 1.74 | 0.88 – 3.42 | 0.109 | 1.72 | 0.91 – 3.23 | 0.094 | |||
| Isolate [C18] | 1.23 | 0.67 – 2.26 | 0.497 | 1.21 | 0.66 – 2.22 | 0.545 | 1.28 | 0.71 – 2.30 | 0.415 | |||
| Isolate [C24] | 1.03 | 0.55 – 1.93 | 0.918 | 0.99 | 0.52 – 1.90 | 0.977 | 1.17 | 0.64 – 2.13 | 0.616 | |||
|
Isolate [C1] × Temperature [linear] |
0.68 | 0.21 – 2.22 | 0.519 | |||||||||
|
Isolate [C18] × Temperature [linear] |
1.70 | 0.59 – 4.91 | 0.324 | |||||||||
|
Isolate [C24] × Temperature [linear] |
0.62 | 0.20 – 1.91 | 0.400 | |||||||||
|
Isolate [C1] × Temperature [quadratic] |
3.40 | 1.08 – 10.67 | 0.036 | |||||||||
|
Isolate [C18] × Temperature [quadratic] |
2.13 | 0.74 – 6.13 | 0.162 | |||||||||
|
Isolate [C24] × Temperature [quadratic] |
0.94 | 0.31 – 2.85 | 0.913 | |||||||||
| Random Effects | ||||||||||||
| σ2 | 3.29 | 3.29 | 3.29 | 3.29 | ||||||||
| τ00 | 0.28 Batch | 0.28 Batch | 0.34 Batch | 0.23 Batch | ||||||||
| ICC | 0.08 | 0.08 | 0.09 | 0.06 | ||||||||
| N | 50 Batch | 50 Batch | 50 Batch | 50 Batch | ||||||||
| Observations | 385 | 385 | 385 | 385 | ||||||||
| Marginal R2 / Conditional R2 | 0.066 / 0.139 | 0.076 / 0.149 | 0.117 / 0.201 | 0.010 / 0.073 | ||||||||
# Step 3: Final model selection based on AICc and interpretability
Inf_final = glmmTMB(data = crushed, InfectedStatus ~ Temperature + (1 | Batch), family = "binomial", REML = TRUE)
# Validating final model with residuals simulation
simulateResiduals(Inf_final, plot = TRUE)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2662472 0.337641 0.2329939 0.1066909 0.2742451 0.0246008 0.5366342 0.9108192 0.1563907 0.003480057 0.5733479 0.4877609 0.14512 0.3991937 0.304557 0.4133086 0.2767652 0.4752524 0.6629793 0.8902504 ...
# Summarizing the final model without intercept
tab_model(Inf_final, show.intercept = FALSE)
| InfectedStatus | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| Temperature [linear] | 0.44 | 0.30 – 0.63 | <0.001 |
| Temperature [quadratic] | 1.03 | 0.71 – 1.49 | 0.872 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 Batch | 0.30 | ||
| ICC | 0.08 | ||
| N Batch | 50 | ||
| Observations | 385 | ||
| Marginal R2 / Conditional R2 | 0.061 / 0.140 | ||
Higher temperatures endured by spores significantly reduce infection
rates, with the odds of infection decreasing by 57% per temperature unit
increase (OR = 0.44, p < 0.001). No intermediate optimal temperature
effect was observed, which suggests a consistent linear decrease.
Although batch effects accounted for some variability, temperature
remains the main driver.
This section explore how temperature and the parasite’s genotype contribute to the survival of infected individuals (also known as virulence) and uninfected individuals (potentially related to defense costs).
surv_obj = Surv(time = infected$DeathAge, event = infected$statusDeath)
fit_cox_5 = coxph(surv_obj ~ Isolate * Temperature, data = infected)
fit_cox_5R = coxme(surv_obj ~ Isolate * Temperature + (1|Batch), data = infected)
AICc(fit_cox_5, fit_cox_5R) # Random effect is significantly improving our results
## df AICc
## fit_cox_5 11.00000 1329.365
## fit_cox_5R 28.43401 1326.367
fit_cox_1 = coxme(surv_obj ~ 1 + (1|Batch), data = infected)
fit_cox_2 = coxme(surv_obj ~ Temperature + (1|Batch), data = infected)
fit_cox_3 = coxme(surv_obj ~ Isolate + (1|Batch), data = infected)
fit_cox_4 = coxme(surv_obj ~ Isolate + Temperature + (1|Batch), data = infected)
fit_cox_5 = coxme(surv_obj ~ Isolate * Temperature + (1|Batch), data = infected)
AICc(fit_cox_1,fit_cox_2,fit_cox_3,fit_cox_4,fit_cox_5)
## df AICc
## fit_cox_1 13.67456 1328.492
## fit_cox_2 16.10514 1331.733
## fit_cox_3 16.71636 1316.222
## fit_cox_4 19.33658 1319.417
## fit_cox_5 28.43401 1326.367
summary(fit_cox_3)
## Mixed effects coxme model
## Formula: surv_obj ~ Isolate + (1 | Batch)
## Data: infected
##
## events, n = 153, 167
##
## Random effects:
## group variable sd variance
## 1 Batch Intercept 0.3965349 0.1572399
## Chisq df p AIC BIC
## Integrated loglik 21.73 4.00 2.270e-04 13.73 1.61
## Penalized loglik 53.57 16.72 9.682e-06 20.14 -30.52
##
## Fixed effects:
## coef exp(coef) se(coef) z p
## IsolateC1 0.5158 1.6751 0.2618 1.97 0.0488
## IsolateC18 -0.5384 0.5837 0.2580 -2.09 0.0369
## IsolateC24 -0.1406 0.8688 0.2588 -0.54 0.5869
### Test assumptions of proportional hazards
(test_ph = cox.zph(fit_cox_3))
## chisq df p
## Isolate 6.15 3 0.1
## GLOBAL 6.15 3 0.1
plot(test_ph)
# Extract results from summary
summary_results = summary(fit_cox_3)
# Extract coefficients, standard errors, and p-values
coef = summary_results$coefficients[, "coef"]
se_coef = summary_results$coefficients[, "se(coef)"]
p_values = summary_results$coefficients[, "p"]
# Calculate HRs and confidence intervals
HR = exp(coef)
lower_CI = exp(coef - 1.96 * se_coef)
upper_CI = exp(coef + 1.96 * se_coef)
# Combine results into a data frame
data.frame(
Predictors = rownames(summary_results$coefficients),
HR = HR,
Lower_95_CI = lower_CI,
Upper_95_CI = upper_CI,
p_value = p_values
)
## Predictors HR Lower_95_CI Upper_95_CI p_value
## IsolateC1 IsolateC1 1.6750547 1.0026986 2.7982569 0.04880475
## IsolateC18 IsolateC18 0.5837033 0.3520170 0.9678783 0.03692993
## IsolateC24 IsolateC24 0.8688436 0.5232149 1.4427900 0.58690328
surv_obj = Surv(time = not_infected$DeathAge, event = not_infected$statusDeath)
fit_cox_5 = coxph(surv_obj ~ Isolate * Temperature, data = not_infected)
fit_cox_5R = coxme(surv_obj ~ Isolate * Temperature + (1|Batch), data = not_infected)
AICc(fit_cox_5, fit_cox_5R) # Random effect
## df AICc
## fit_cox_5 11.00000 1125.461
## fit_cox_5R 27.43641 1130.994
fit_cox_1 = coxph(surv_obj ~ 1 , data = not_infected)
fit_cox_2 = coxph(surv_obj ~ Temperature, data = not_infected)
fit_cox_3 = coxph(surv_obj ~ Isolate , data = not_infected)
fit_cox_4 = coxph(surv_obj ~ Isolate + Temperature , data = not_infected)
fit_cox_5 = coxph(surv_obj ~ Isolate * Temperature , data = not_infected)
AICc(fit_cox_1,fit_cox_2,fit_cox_3,fit_cox_4,fit_cox_5)
## df AICc
## fit_cox_1 0 1117.001
## fit_cox_2 2 1113.668
## fit_cox_3 3 1122.686
## fit_cox_4 5 1119.611
## fit_cox_5 11 1125.461
summary(fit_cox_2)
## Call:
## coxph(formula = surv_obj ~ Temperature, data = not_infected)
##
## n= 218, number of events= 110
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Temperature.L -0.4530 0.6357 0.1633 -2.773 0.00555 **
## Temperature.Q 0.1291 1.1379 0.1705 0.757 0.44882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Temperature.L 0.6357 1.5730 0.4616 0.8756
## Temperature.Q 1.1379 0.8789 0.8146 1.5894
##
## Concordance= 0.577 (se = 0.026 )
## Likelihood ratio test= 7.45 on 2 df, p=0.02
## Wald test = 7.98 on 2 df, p=0.02
## Score (logrank) test = 8.21 on 2 df, p=0.02
anova(fit_cox_1, fit_cox_2)
## Analysis of Deviance Table
## Cox model: response is surv_obj
## Model 1: ~ 1
## Model 2: ~ Temperature
## loglik Chisq Df Pr(>|Chi|)
## 1 -558.50
## 2 -554.78 7.4457 2 0.02416 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Test assumptions of proportional hazards
(test_ph = cox.zph(fit_cox_2))
## chisq df p
## Temperature 6.48 2 0.039
## GLOBAL 6.48 2 0.039
plot(test_ph) ### Not so bad
### Weibull approach to compare results
fit_weibull1 = survreg(surv_obj ~ 1 , data = not_infected, dist = "weibull")
fit_weibull2 = survreg(surv_obj ~ Temperature, data = not_infected, dist = "weibull")
fit_weibull3 = survreg(surv_obj ~ Isolate, data = not_infected, dist = "weibull")
fit_weibull4 = survreg(surv_obj ~ Isolate + Temperature, data = not_infected, dist = "weibull")
fit_weibull5 = survreg(surv_obj ~ Isolate * Temperature, data = not_infected, dist = "weibull")
AICc(fit_weibull1,fit_weibull2,fit_weibull3,fit_weibull4,fit_weibull5)
## df AICc
## fit_weibull1 2 1212.177
## fit_weibull2 4 1207.967
## fit_weibull3 5 1217.837
## fit_weibull4 7 1213.874
## fit_weibull5 13 1219.663
# Standardi@zed residuals plot
residuals_std = residuals(fit_weibull2, type = "matrix")
qqnorm(residuals_std)
qqline(residuals_std, col = "red")
summary(fit_weibull2) # Same results, temperature is the main effect, with a significant decrease in mortality risk at higher temperatures (p = 0.003).
##
## Call:
## survreg(formula = surv_obj ~ Temperature, data = not_infected,
## dist = "weibull")
## Value Std. Error z p
## (Intercept) 4.3576 0.0589 74.01 < 2e-16
## Temperature.L 0.2701 0.0931 2.90 0.0037
## Temperature.Q -0.0777 0.0960 -0.81 0.4182
## Log(scale) -0.5755 0.0864 -6.66 2.7e-11
##
## Scale= 0.562
##
## Weibull distribution
## Loglik(model)= -599.9 Loglik(intercept only)= -604.1
## Chisq= 8.34 on 2 degrees of freedom, p= 0.015
## Number of Newton-Raphson Iterations: 5
## n= 218
# So I decided to keep the cox model despite a slight violation of the proportional hazard assumption, as it is more interpretable and the results are consistent with the Weibull model.
# Extract results from summary
summary_results = summary(fit_cox_2)
# Extract coefficients, standard errors, and p-values
coef = summary_results$coefficients[, "coef"]
se_coef = summary_results$coefficients[, "se(coef)"]
p_values = summary_results$coefficients[, "Pr(>|z|)"]
# Calculate HRs and confidence intervals
HR = exp(coef)
lower_CI = exp(coef - 1.96 * se_coef)
upper_CI = exp(coef + 1.96 * se_coef)
# Combine results into a data frame
data.frame(
Predictors = rownames(summary_results$coefficients),
HR = HR,
Lower_95_CI = lower_CI,
Upper_95_CI = upper_CI,
p_value = p_values
)
## Predictors HR Lower_95_CI Upper_95_CI p_value
## Temperature.L Temperature.L 0.6357221 0.4615649 0.8755921 0.005547673
## Temperature.Q Temperature.Q 1.1378501 0.8146002 1.5893723 0.448821505
The model indicates that temperature significantly decreases the hazard
rate, with a 37.1% reduction in mortality risk (HR = 0.63, p = 0.006),
suggesting potential temperature-mediated effects on host survival.
## C1 C14 C18 C24 CTRL
## 45 58 67 60 27
We start first by running an analysis including the husbandry control along with the other genotype.
Warning! It is important to note that under these conditions it is not possible to test the interaction between heatwave temperature and the absence of spores, as our husbandry controls were infected with crushed Daphnia bodies that did not experience any heatwaves.
surv_obj = Surv(time = dataC$DeathAge, event = dataC$statusDeath)
fit_cox_5 = coxph(surv_obj ~ Temperature + Isolate, data = dataC)
fit_cox_5R = coxme(surv_obj ~ Temperature + Isolate + (1|Batch), data = dataC)
# Cannot evaluate the interaction as the unexposed individuals do not have a heatwave condition (no spores)
AICc(fit_cox_5, fit_cox_5R) # No random effect
## df AICc
## fit_cox_5 6.00000 1376.118
## fit_cox_5R 12.21118 1376.970
fit_cox_1 = coxph(surv_obj ~ 1, data = dataC)
fit_cox_2 = coxph(surv_obj ~ Temperature, data = dataC)
fit_cox_3 = coxph(surv_obj ~ Isolate, data = dataC)
fit_cox_4 = coxph(surv_obj ~ Isolate + Temperature, data = dataC)
AICc(fit_cox_1,fit_cox_2,fit_cox_3,fit_cox_4)
## df AICc
## fit_cox_1 0 1376.198
## fit_cox_2 2 1368.589
## fit_cox_3 4 1379.257
## fit_cox_4 6 1376.118
### Test assumptions of proportional hazards
test_ph = cox.zph(fit_cox_2)
test_ph
## chisq df p
## Temperature 4.38 2 0.11
## GLOBAL 4.38 2 0.11
plot(test_ph) # Not the best. fit
### Weibull approach to compare results
fit_weibull1 = survreg(surv_obj ~ 1 , data = dataC, dist = "weibull")
fit_weibull2 = survreg(surv_obj ~ Temperature, data = dataC, dist = "weibull")
fit_weibull3 = survreg(surv_obj ~ Isolate, data = dataC, dist = "weibull")
fit_weibull4 = survreg(surv_obj ~ Isolate + Temperature, data = dataC, dist = "weibull")
fit_weibull5 = survreg(surv_obj ~ Isolate * Temperature, data = dataC, dist = "weibull")
AICc(fit_weibull1,fit_weibull2,fit_weibull3,fit_weibull4,fit_weibull5)
## df AICc
## fit_weibull1 2 1438.865
## fit_weibull2 4 1429.968
## fit_weibull3 6 1441.145
## fit_weibull4 8 1437.353
## fit_weibull5 16 1447.669
# Standardi@zed residuals plot
residuals_std = residuals(fit_weibull2, type = "matrix")
qqnorm(residuals_std)
qqline(residuals_std, col = "red")
summary(fit_weibull2) # not much much better, but sit's also because the factor are not equilibrated as the design is not meant to test those hypothesese.
##
## Call:
## survreg(formula = surv_obj ~ Temperature, data = dataC, dist = "weibull")
## Value Std. Error z p
## (Intercept) 4.3525 0.0523 83.27 < 2e-16
## Temperature.L 0.2662 0.0783 3.40 0.00067
## Temperature.Q -0.1057 0.0888 -1.19 0.23395
## Log(scale) -0.6169 0.0782 -7.89 3.1e-15
##
## Scale= 0.54
##
## Weibull distribution
## Loglik(model)= -710.9 Loglik(intercept only)= -717.4
## Chisq= 13.01 on 2 degrees of freedom, p= 0.0015
## Number of Newton-Raphson Iterations: 5
## n= 257
summary(fit_cox_2)
## Call:
## coxph(formula = surv_obj ~ Temperature, data = dataC)
##
## n= 257, number of events= 132
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Temperature.L -0.4683 0.6261 0.1435 -3.264 0.0011 **
## Temperature.Q 0.1847 1.2028 0.1644 1.123 0.2614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Temperature.L 0.6261 1.5973 0.4726 0.8294
## Temperature.Q 1.2028 0.8314 0.8715 1.6601
##
## Concordance= 0.582 (se = 0.023 )
## Likelihood ratio test= 11.7 on 2 df, p=0.003
## Wald test = 12.24 on 2 df, p=0.002
## Score (logrank) test = 12.64 on 2 df, p=0.002
Secondly, we analysed the survival of uninfected Daphnias only at 20°C, as this is a condition were the spores endured no heatwave, like our husbandry controls.
## df AICc
## fit_cox_5 4.00000 433.5022
## fit_cox_5R 4.01947 433.5129
## df AICc
## fit_cox_1 0 429.8015
## fit_cox_2 4 433.5022
## Call:
## coxph(formula = surv_obj ~ Isolate, data = dataC20)
##
## n= 84, number of events= 54
##
## coef exp(coef) se(coef) z Pr(>|z|)
## IsolateC1 0.7564 2.1306 0.4438 1.704 0.0883 .
## IsolateC14 -0.1502 0.8605 0.4052 -0.371 0.7108
## IsolateC18 -0.2094 0.8111 0.3910 -0.536 0.5923
## IsolateC24 -0.4021 0.6689 0.4058 -0.991 0.3218
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## IsolateC1 2.1306 0.4694 0.8928 5.084
## IsolateC14 0.8605 1.1621 0.3889 1.904
## IsolateC18 0.8111 1.2329 0.3769 1.745
## IsolateC24 0.6689 1.4949 0.3020 1.482
##
## Concordance= 0.571 (se = 0.044 )
## Likelihood ratio test= 5.12 on 4 df, p=0.3
## Wald test = 5.98 on 4 df, p=0.2
## Score (logrank) test = 6.39 on 4 df, p=0.2
## chisq df p
## Isolate 4.19 4 0.38
## GLOBAL 4.19 4 0.38
We can see that there is no significant difference between the model
comparing the survival of each spore lineage (including the group
without spores) and the null model. So it seem than at 20°C, in stable
conditions, the exposure to dead Daphnias or to spores from different
genotype influence similarly the survival of non infected
individuals.
In conclusion, when we included the husbandry control in our analysis, comparing the additive effects of heatwave temperature and different spores genotypes (including non spores at all), the main effect explaining survival variation is still temperature.
In this section we explored how the exposure and the infection by heat-stressed parasites affected the host reproduction. First we evaluated the number of individuals fully castrated in exposed uninfected and in infected Daphnias. This divisions into two sub analysis helped us to find more adapted models to describe the host reproductive patern.
One of the consequences of Pasteuria ramosa infection in Daphnia is the reduction in allocation to reproduction, which can also potentially be due to the energetic cost of immune response for uninfected hosts. Here we focused on the proportion of non-reproducing individuals.
FC1 = glmmTMB(data=infected, Never_repro ~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=infected, Never_repro ~ Isolate + Temperature + DeathAge + (1|Batch), family = "binomial", REML=T)
AICc(FC1, FC1R)
## df AICc
## FC1 7 239.8814
## FC1R 8 238.7983
FC1 = glmmTMB(data=infected, Never_repro ~ Isolate * Temperature + DeathAge, family = "binomial",
na.action = na.fail)
model_selection = dredge(FC1)
# see best models
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge,
## data = infected, family = "binomial", na.action = na.fail,
## ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 16 0.79900 + -0.02228 + + + 13 -98.284
## 15 -0.19180 + + + + 12 -99.622
## 4 0.98420 + -0.01982 + 5 -108.183
## 3 0.11780 + + 4 -109.410
## 2 0.97250 + -0.02725 2 -111.825
## 8 0.94670 + -0.01938 + + 7 -107.700
## 7 0.09487 + + + 6 -108.859
## 1 -0.25280 + 1 -114.432
## 6 0.91880 + -0.02678 + 4 -111.381
## 5 -0.28880 + + 3 -113.875
## AICc delta weight
## 16 224.9 0.00 0.321
## 15 225.3 0.32 0.273
## 4 226.7 1.79 0.131
## 3 227.1 2.12 0.111
## 2 227.7 2.78 0.080
## 8 230.1 5.16 0.024
## 7 230.2 5.30 0.023
## 1 230.9 5.94 0.016
## 6 231.0 6.06 0.016
## 5 233.9 8.95 0.004
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| Never_repro | Never_repro | Never_repro | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| DeathAge | 0.98 | 0.95 – 1.00 | 0.104 | 0.98 | 0.96 – 1.01 | 0.120 | |||
| Isolate [C1] | 2.27 | 0.66 – 7.83 | 0.193 | 2.46 | 0.72 – 8.38 | 0.151 | 1.23 | 0.49 – 3.10 | 0.658 |
| Isolate [C18] | 0.55 | 0.18 – 1.69 | 0.298 | 0.50 | 0.16 – 1.49 | 0.212 | 0.44 | 0.18 – 1.11 | 0.083 |
| Isolate [C24] | 0.65 | 0.20 – 2.09 | 0.472 | 0.66 | 0.21 – 2.08 | 0.474 | 0.50 | 0.20 – 1.28 | 0.150 |
| Temperature [linear] | 0.17 | 0.03 – 0.96 | 0.044 | 0.14 | 0.03 – 0.77 | 0.024 | |||
| Temperature [quadratic] | 0.28 | 0.07 – 1.12 | 0.073 | 0.29 | 0.07 – 1.16 | 0.080 | |||
|
Isolate [C1] × Temperature [linear] |
32.02 | 3.12 – 328.91 | 0.004 | 36.19 | 3.56 – 367.88 | 0.002 | |||
|
Isolate [C18] × Temperature [linear] |
2.04 | 0.26 – 16.14 | 0.497 | 2.80 | 0.37 – 21.15 | 0.319 | |||
|
Isolate [C24] × Temperature [linear] |
5.08 | 0.56 – 45.81 | 0.147 | 6.46 | 0.74 – 56.47 | 0.092 | |||
|
Isolate [C1] × Temperature [quadratic] |
10.98 | 1.60 – 75.39 | 0.015 | 10.14 | 1.50 – 68.66 | 0.018 | |||
|
Isolate [C18] × Temperature [quadratic] |
3.86 | 0.64 – 23.40 | 0.142 | 3.51 | 0.59 – 20.92 | 0.168 | |||
|
Isolate [C24] × Temperature [quadratic] |
3.68 | 0.59 – 23.04 | 0.163 | 3.31 | 0.54 – 20.28 | 0.196 | |||
| Observations | 167 | 167 | 167 | ||||||
| R2 Tjur | 0.178 | 0.163 | 0.074 | ||||||
best_model = get.models(model_selection, 2)[[1]]
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3977238 0.8536237 0.8579498 0.6620259 0.61213 0.9357635 0.05493602 0.2030844 0.1008657 0.9558345 0.1985479 0.8851658 0.4544661 0.1936561 0.1833452 0.2271874 0.5898748 0.8261575 0.5827213 0.7709696 ...
tab_model(best_model, show.intercept = F)
| Never_repro | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| Isolate [C1] | 2.46 | 0.72 – 8.38 | 0.151 |
| Isolate [C18] | 0.50 | 0.16 – 1.49 | 0.212 |
| Isolate [C24] | 0.66 | 0.21 – 2.08 | 0.474 |
| Temperature [linear] | 0.14 | 0.03 – 0.77 | 0.024 |
| Temperature [quadratic] | 0.29 | 0.07 – 1.16 | 0.080 |
|
Isolate [C1] × Temperature [linear] |
36.19 | 3.56 – 367.88 | 0.002 |
|
Isolate [C18] × Temperature [linear] |
2.80 | 0.37 – 21.15 | 0.319 |
|
Isolate [C24] × Temperature [linear] |
6.46 | 0.74 – 56.47 | 0.092 |
|
Isolate [C1] × Temperature [quadratic] |
10.14 | 1.50 – 68.66 | 0.018 |
|
Isolate [C18] × Temperature [quadratic] |
3.51 | 0.59 – 20.92 | 0.168 |
|
Isolate [C24] × Temperature [quadratic] |
3.31 | 0.54 – 20.28 | 0.196 |
| Observations | 167 | ||
| R2 Tjur | 0.163 | ||
## `summarise()` has grouped output by 'Temperature'. You can override using the
## `.groups` argument.
A heatwave experienced by the parasite appears to relax its castrating
effect on the host for most genotypes, with a strong opposit effect for
C1 genotyp, with a lower probability of total sterility in genotype C1,
which shows a dramatically increased risk of host sterility after
exposure to 40°C. This suggests an interaction between genotype and
heatwave experienced by the parasite.
options(na.action = "na.fail")
FC1 = glmmTMB(data=not_infected, Never_repro~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=not_infected, Never_repro ~ Isolate + Temperature+ DeathAge + (1|Batch), family = "binomial", REML=T)
AICc(FC1, FC1R)
## df AICc
## FC1 7 118.8705
## FC1R 8 106.9146
FC1 = glmmTMB(data=not_infected, Never_repro ~ Isolate * Temperature + DeathAge + (1|Batch), family = "binomial")
model_selection = dredge(FC1)
# Table of model comparison
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge +
## (1 | Batch), data = not_infected, family = "binomial", ziformula = ~0,
## dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 2 1.672 + -0.1281 3 -48.375
## 6 1.459 + -0.1208 + 5 -46.911
## 4 2.324 + -0.1529 + 6 -46.318
## 8 2.030 + -0.1347 + + 8 -45.341
## 16 13.600 + -0.8717 + + + 14 -40.057
## 5 -2.909 + + 4 -67.021
## 7 -2.961 + + + 7 -65.982
## 1 -2.830 + 2 -72.618
## 15 -3.165 + + + + 13 -61.318
## 3 -2.927 + + 5 -71.734
## AICc delta weight
## 2 102.9 0.00 0.499
## 6 104.1 1.24 0.268
## 4 105.0 2.17 0.168
## 8 107.4 4.51 0.052
## 16 110.2 7.32 0.013
## 5 142.2 39.37 0.000
## 7 146.5 43.63 0.000
## 1 149.3 46.43 0.000
## 15 150.4 47.56 0.000
## 3 153.8 50.89 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | Batch)
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3977238 0.8536237 0.8579498 0.6620259 0.61213 0.9357635 0.05493602 0.2030844 0.1008657 0.9558345 0.1985479 0.8851658 0.4544661 0.1936561 0.1833452 0.2271874 0.5898748 0.8261575 0.5827213 0.7709696 ...
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| Never_repro | Never_repro | |||||
|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p |
| DeathAge | 0.88 | 0.82 – 0.94 | <0.001 | 0.89 | 0.83 – 0.95 | 0.001 |
| Temperature [linear] | 0.36 | 0.10 – 1.30 | 0.120 | |||
| Temperature [quadratic] | 0.61 | 0.17 – 2.17 | 0.442 | |||
| Random Effects | ||||||
| σ2 | 3.29 | 3.29 | ||||
| τ00 | 6.14 Batch | 5.52 Batch | ||||
| ICC | 0.65 | 0.63 | ||||
| N | 50 Batch | 50 Batch | ||||
| Observations | 218 | 218 | ||||
| Marginal R2 / Conditional R2 | 0.385 / 0.786 | 0.421 / 0.784 | ||||
# best model
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
## Family: binomial ( logit )
## Formula: Never_repro ~ DeathAge + (1 | Batch)
## Data: not_infected
##
## AIC BIC logLik deviance df.resid
## 102.8 112.9 -48.4 96.8 215
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Batch (Intercept) 6.137 2.477
## Number of obs: 218, groups: Batch, 50
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.67209 0.98865 1.691 0.09078 .
## DeathAge -0.12808 0.03582 -3.575 0.00035 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(best_model)
| Never_repro | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 5.32 | 0.77 – 36.96 | 0.091 |
| DeathAge | 0.88 | 0.82 – 0.94 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 Batch | 6.14 | ||
| ICC | 0.65 | ||
| N Batch | 50 | ||
| Observations | 218 | ||
| Marginal R2 / Conditional R2 | 0.385 / 0.786 | ||
Only the age at death seems to reprduce at least once in their life
time.
## `summarise()` has grouped output by 'DeathAgeGroup'. You can override using the
## `.groups` argument.
options(na.action = "na.fail")
dataC_ninf = subset(dataC, InfectedStatus!="infected")
summary(dataC_ninf)
## ID Temperature Isolate Clutch Nb_offspring
## Min. : 2 20: 84 C1 :45 Min. : 1.000 Min. : 0.00
## 1st Qu.: 218 30: 73 C14 :58 1st Qu.: 3.750 1st Qu.: 10.00
## Median : 668 40:100 C18 :67 Median : 8.000 Median : 47.00
## Mean : 588 C24 :60 Mean : 7.368 Mean : 43.02
## 3rd Qu.: 969 CTRL:27 3rd Qu.:11.000 3rd Qu.: 68.00
## Max. :1100 Max. :16.000 Max. :111.00
## NA's :29
## FirstClutchAge LastClutchAge DeathAge Accidental_death
## Min. : 1.00 Min. :10.00 Min. :15.00 Min. :1
## 1st Qu.:14.00 1st Qu.:39.00 1st Qu.:31.00 1st Qu.:1
## Median :18.00 Median :59.00 Median :62.00 Median :1
## Mean :21.77 Mean :50.47 Mean :51.47 Mean :1
## 3rd Qu.:27.00 3rd Qu.:67.00 3rd Qu.:67.00 3rd Qu.:1
## Max. :67.00 Max. :67.00 Max. :67.00 Max. :1
## NA's :30 NA's :32 NA's :250
## Volume_counted InfectedStatus Number of mature spores
## Min. :150 Length:257 Length:257
## 1st Qu.:150 Class :character Class :character
## Median :150 Mode :character Mode :character
## Mean :152
## 3rd Qu.:150
## Max. :200
## NA's :12
## Number of imature spores Nb of white spores cauliflowers grapes
## Length:257 Length:257 no :236 no :236
## Class :character Class :character NA's: 21 NA's: 21
## Mode :character Mode :character
##
##
##
##
## whites group statusDeath Number of immature spores
## Min. :0 C1440 : 29 Min. :0.0000 Length:257
## 1st Qu.:0 C2440 : 28 1st Qu.:0.0000 Class :character
## Median :0 NACTRL : 27 Median :1.0000 Mode :character
## Mean :0 C1840 : 26 Mean :0.5136
## 3rd Qu.:0 C1830 : 24 3rd Qu.:1.0000
## Max. :0 C130 : 19 Max. :1.0000
## (Other):104
## Never_repro Batch imm mat mature_load
## Min. :0.0000 2 : 8 Min. :0 Min. :0 Min. :0
## 1st Qu.:0.0000 17 : 8 1st Qu.:0 1st Qu.:0 1st Qu.:0
## Median :0.0000 28 : 8 Median :0 Median :0 Median :0
## Mean :0.1128 31 : 8 Mean :0 Mean :0 Mean :0
## 3rd Qu.:0.0000 3 : 7 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
## Max. :1.0000 6 : 7 Max. :0 Max. :0 Max. :0
## (Other):211 NA's :12
## imm_load ratio repro_rate stade_max
## Min. :0 Min. : NA Min. :0.0000 None:257
## 1st Qu.:0 1st Qu.: NA 1st Qu.:0.3095
## Median :0 Median : NA Median :0.8358
## Mean :0 Mean :NaN Mean :0.7206
## 3rd Qu.:0 3rd Qu.: NA 3rd Qu.:1.0597
## Max. :0 Max. : NA Max. :1.6567
## NA's :12 NA's :257
FC1 = glmmTMB(data=dataC_ninf, Never_repro~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=dataC_ninf, Never_repro ~ Isolate + Temperature+ DeathAge + (1|Batch), family = "binomial", REML=T)
AICc(FC1, FC1R)
## df AICc
## FC1 8 134.3687
## FC1R 9 125.3196
FC1 = glmmTMB(data=not_infected, Never_repro ~ Isolate * Temperature + DeathAge + (1|Batch), family = "binomial")
model_selection = dredge(FC1)
# Table of model comparison
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge +
## (1 | Batch), data = not_infected, family = "binomial", ziformula = ~0,
## dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 2 1.672 + -0.1281 3 -48.375
## 6 1.459 + -0.1208 + 5 -46.911
## 4 2.324 + -0.1529 + 6 -46.318
## 8 2.030 + -0.1347 + + 8 -45.341
## 16 13.600 + -0.8717 + + + 14 -40.057
## 5 -2.909 + + 4 -67.021
## 7 -2.961 + + + 7 -65.982
## 1 -2.830 + 2 -72.618
## 15 -3.165 + + + + 13 -61.318
## 3 -2.927 + + 5 -71.734
## AICc delta weight
## 2 102.9 0.00 0.499
## 6 104.1 1.24 0.268
## 4 105.0 2.17 0.168
## 8 107.4 4.51 0.052
## 16 110.2 7.32 0.013
## 5 142.2 39.37 0.000
## 7 146.5 43.63 0.000
## 1 149.3 46.43 0.000
## 15 150.4 47.56 0.000
## 3 153.8 50.89 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | Batch)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| Never_repro | Never_repro | |||||
|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p |
| DeathAge | 0.88 | 0.82 – 0.94 | <0.001 | 0.89 | 0.83 – 0.95 | 0.001 |
| Temperature [linear] | 0.36 | 0.10 – 1.30 | 0.120 | |||
| Temperature [quadratic] | 0.61 | 0.17 – 2.17 | 0.442 | |||
| Random Effects | ||||||
| σ2 | 3.29 | 3.29 | ||||
| τ00 | 6.14 Batch | 5.52 Batch | ||||
| ICC | 0.65 | 0.63 | ||||
| N | 50 Batch | 50 Batch | ||||
| Observations | 218 | 218 | ||||
| Marginal R2 / Conditional R2 | 0.385 / 0.786 | 0.421 / 0.784 | ||||
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
## Family: binomial ( logit )
## Formula: Never_repro ~ DeathAge + (1 | Batch)
## Data: not_infected
##
## AIC BIC logLik deviance df.resid
## 102.8 112.9 -48.4 96.8 215
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Batch (Intercept) 6.137 2.477
## Number of obs: 218, groups: Batch, 50
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.67209 0.98865 1.691 0.09078 .
## DeathAge -0.12808 0.03582 -3.575 0.00035 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(best_model)
| Never_repro | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 5.32 | 0.77 – 36.96 | 0.091 |
| DeathAge | 0.88 | 0.82 – 0.94 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 Batch | 6.14 | ||
| ICC | 0.65 | ||
| N Batch | 50 | ||
| Observations | 218 | ||
| Marginal R2 / Conditional R2 | 0.385 / 0.786 | ||
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.9088821 0.9008778 0.9317442 0.7483929 0.241452 0.2361889 0.7475164 0.4387352 0.390151 0.5425091 0.609537 0.7492443 0.6104253 0.7648734 0.5350931 0.2879339 0.134203 0.1612245 0.3157523 0.1570396 ...
Only the age at death seems to reprduce at least once in their life
time, the husbandry control inclusion does not change this.
## `summarise()` has grouped output by 'DeathAgeGroup'. You can override using the
## `.groups` argument.
In this section we were interested in understanding the role of parasite heat-stress and genotype on the number of offspring produced. Indeed, one of the consequences of infection in Daphnia is the reduction in allocation to parasite reproduction, (or a potential cost of immune response for uninfected individuals), even if the host is not fully castrated.
repro_inf = subset(data, data$Clutch>0 & InfectedStatus=="infected" & Temperature!="NA")
repro_inf$Temperature = factor(repro_inf$Temperature, ordered = T)
FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "nbinom2", REML=T)
#FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "poisson", REML=T)
#FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "gaussian", REML=T)
simulateResiduals(FC1, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.8741857 0.7996104 0.4441007 0.4213267 0.3856629 0.2460142 0.4371622 0.9950848 0.5756341 0.3335869 0.1840254 0.6387504 0.2376711 0.8736609 0.9228149 0.7387757 0.2361033 0.4922073 0.4105355 0.1956212 ...
FC1R = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge + (1|Batch), family = "nbinom2", REML=T)
AICc(FC1, FC1R)# no batch effect
## df AICc
## FC1 8 601.8980
## FC1R 9 601.5721
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge ,
data = repro_inf, family = "nbinom2", na.action = na.fail)
model_selection = dredge(FC1)
model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge,
## data = repro_inf, family = "nbinom2", na.action = na.fail,
## ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 4 0.9519 + 0.03153 + 6 -283.958
## 2 0.2723 + 0.03680 3 -289.038
## 8 0.9455 + 0.03176 + + 8 -283.902
## 16 0.7038 + 0.02896 + + + 14 -277.332
## 6 0.3050 + 0.03628 + 5 -288.741
## 3 2.7120 + + 5 -291.710
## 7 2.7390 + + + 7 -291.682
## 15 2.2290 + + + + 13 -284.113
## 1 2.1350 + 2 -300.042
## 5 2.1400 + + 4 -299.296
## AICc delta weight
## 4 580.9 0.00 0.750
## 2 584.3 3.46 0.133
## 8 585.5 4.62 0.075
## 16 588.0 7.10 0.022
## 6 588.2 7.28 0.020
## 3 594.1 13.22 0.001
## 7 598.7 17.79 0.000
## 15 598.8 17.89 0.000
## 1 604.2 23.34 0.000
## 5 607.0 26.16 0.000
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| Nb_offspring | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| DeathAge | 1.03 | 1.02 – 1.05 | <0.001 |
| Isolate [C1] | 0.44 | 0.22 – 0.89 | 0.021 |
| Isolate [C18] | 0.48 | 0.28 – 0.82 | 0.008 |
| Isolate [C24] | 0.80 | 0.45 – 1.39 | 0.424 |
| Observations | 94 | ||
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.8445634 0.7757172 0.3188654 0.4678265 0.4046863 0.2574613 0.3273075 0.988557 0.6029993 0.2770807 0.2222652 0.5937299 0.200338 0.8887268 0.8925955 0.6897923 0.2627872 0.457771 0.3542096 0.1563993 ...
tab_model(best_model, show.intercept = F)
| Nb_offspring | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| DeathAge | 1.03 | 1.02 – 1.05 | <0.001 |
| Isolate [C1] | 0.44 | 0.22 – 0.89 | 0.021 |
| Isolate [C18] | 0.48 | 0.28 – 0.82 | 0.008 |
| Isolate [C24] | 0.80 | 0.45 – 1.39 | 0.424 |
| Observations | 94 | ||
summary(best_model)
## Family: nbinom2 ( log )
## Formula: Nb_offspring ~ DeathAge + Isolate + 1
## Data: repro_inf
##
## AIC BIC logLik deviance df.resid
## 579.9 595.2 -284.0 567.9 88
##
##
## Dispersion parameter for nbinom2 family (): 1.4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.951947 0.468983 2.030 0.04238 *
## DeathAge 0.031525 0.007622 4.136 3.54e-05 ***
## IsolateC1 -0.812114 0.352549 -2.304 0.02125 *
## IsolateC18 -0.738421 0.276490 -2.671 0.00757 **
## IsolateC24 -0.228773 0.286424 -0.799 0.42445
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
repro_nf = subset(data, data$Clutch>0 & InfectedStatus=="not_inf" & Temperature!="NA")
repro_nf$Temperature = factor(repro_nf$Temperature, ordered = T)
FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature + DeathAge+ (1|Batch), family = "gaussian", REML=T)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature+ DeathAge + (1|Batch), family = "poisson")
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature + DeathAge+ (1|Batch), family = "nbinom2", REML=T)
#FC1 = glmmTMB(data=repro_nf, log(Nb_offspring) ~ Isolate * Temperature + DeathAge + (1|Batch), family = "gaussian", REML=T)
FC1R = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature + DeathAge + (1|Batch), family = "gaussian", REML=F)
FC1 = glmmTMB(Nb_offspring ~ Isolate + Temperature + DeathAge, data = repro_nf, family = "gaussian", REML=F)
AICc(FC1, FC1R)# Batch effect NOT important! But significantly improved the residuals
## df AICc
## FC1 8 1659.382
## FC1R 9 1660.230
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge + (1|Batch), data = repro_nf, family = "gaussian", na.action = na.fail)
model_selection = dredge(FC1)
model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge +
## (1 | Batch), data = repro_nf, family = "gaussian", na.action = na.fail,
## ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 2 -22.63 + 1.306 4 -823.741
## 4 -20.37 + 1.313 + 7 -821.201
## 6 -22.45 + 1.294 + 6 -823.160
## 8 -20.41 + 1.303 + + 9 -820.623
## 5 48.57 + + 5 -913.919
## 1 50.26 + 3 -916.804
## 7 51.89 + + + 8 -913.084
## 3 53.77 + + 6 -915.627
## 15 52.47 + + + + 14 -912.123
## 16 -20.73 + 1.307 + + + 15
## AICc delta
## 2 1655.7 0.00
## 4 1657.0 1.31
## 6 1658.8 3.08
## 8 1660.2 4.54
## 5 1838.2 182.46
## 1 1839.7 184.04
## 7 1843.0 187.26
## 3 1843.7 188.01
## 15 1854.6 198.91
## 16
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | Batch)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
selected_models
## $`2`
## Formula: Nb_offspring ~ DeathAge + (1 | Batch)
## Data: repro_nf
## AIC BIC logLik df.resid
## 1655.4818 1668.5326 -823.7409 189
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## Batch (Intercept) 4.599
## Residual 16.721
##
## Number of obs: 193 / Conditional model: Batch, 50
##
## Dispersion estimate for gaussian family (sigma^2): 280
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) DeathAge
## -22.631 1.306
##
## $`4`
## Formula: Nb_offspring ~ DeathAge + Isolate + (1 | Batch)
## Data: repro_nf
## AIC BIC logLik df.resid
## 1656.4026 1679.2414 -821.2013 186
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## Batch (Intercept) 4.633
## Residual 16.482
##
## Number of obs: 193 / Conditional model: Batch, 50
##
## Dispersion estimate for gaussian family (sigma^2): 272
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) DeathAge IsolateC1 IsolateC18 IsolateC24
## -20.3723 1.3133 -4.2733 -0.6133 -6.5511
##
## $<NA>
## NULL
##
## attr(,"rank")
## function (x)
## do.call("rank", list(x))
## <environment: 0x137323810>
## attr(,"call")
## AICc(x)
## attr(,"class")
## [1] "function" "rankFunction"
## attr(,"beta")
## [1] "none"
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.78 0.656 0.584 0.568 0.468 0.564 0.596 0.448 0.512 0.388 0.464 0.368 0.504 0.452 0.672 0.536 0.708 0.592 0.356 0.488 ...
tab_model(best_model, show.intercept = F) # No effect
| Nb_offspring | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| DeathAge | 1.31 | 1.16 – 1.45 | <0.001 |
| Random Effects | |||
| σ2 | 279.59 | ||
| τ00 Batch | 21.15 | ||
| ICC | 0.07 | ||
| N Batch | 50 | ||
| Observations | 193 | ||
| Marginal R2 / Conditional R2 | 0.622 / 0.649 | ||
summary(best_model)
## Family: gaussian ( identity )
## Formula: Nb_offspring ~ DeathAge + (1 | Batch)
## Data: repro_nf
##
## AIC BIC logLik deviance df.resid
## 1655.5 1668.5 -823.7 1647.5 189
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Batch (Intercept) 21.15 4.599
## Residual 279.59 16.721
## Number of obs: 193, groups: Batch, 50
##
## Dispersion estimate for gaussian family (sigma^2): 280
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.63107 4.34365 -5.21 1.89e-07 ***
## DeathAge 1.30586 0.07379 17.70 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculation of daily reproduction rate
repro_nf$daily_repro_rate = repro_nf$Nb_offspring / repro_nf$DeathAge
mean(repro_nf$daily_repro_rate, na.rm = TRUE)
## [1] 0.8219201
(sem = sd(repro_nf$daily_repro_rate, na.rm = TRUE) / sqrt(length(repro_nf$daily_repro_rate)))
## [1] 0.02774078
Only the effect of time is important here; each day, a Daphnia produces on average 0.30 more offspring.
dataC$Isolate = as.character(dataC$Isolate)
dataC$Isolate[is.na(dataC$Isolate)] = "unexposed"
dataC$Isolate = as.factor(dataC$Isolate)
dataC$Temperature = as.character(dataC$Temperature)
dataC$Temperature[dataC$Temperature == "CTRL"] = 20
dataC$Temperature = factor(dataC$Temperature, ordered = TRUE)
repro_nfC = subset(dataC, dataC$Clutch>0 & InfectedStatus!="inf")
repro_nfC$Temperature = factor(repro_nfC$Temperature, labels = c("20", "30","40"), ordered = T)
FC1 = glmmTMB(data=repro_nfC, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "gaussian", REML=F)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature+ DeathAge, family = "poisson", REML=F)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "nbinom2", REML=T)
simulateResiduals(FC1, plot=T) # I don't find better and I have no reason to exclude the outlier.
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.724 0.608 0.58 0.488 0.6 0.48 0.656 0.54 0.672 0.492 0.484 0.408 0.4 0.548 0.464 0.352 0.8 0.624 0.824 0.668 ...
FC1R = glmmTMB(data=repro_nfC, Nb_offspring ~ Isolate + Temperature + DeathAge + (1|Batch), family = "gaussian", REML=F)
FC1 = glmmTMB(Nb_offspring ~ Isolate + Temperature + DeathAge, data = repro_nfC, family = "gaussian", REML=F)
AICc(FC1, FC1R)# Batch effect NOT important and do not improve residuals anyway.
## df AICc
## FC1 9 1945.748
## FC1R 10 1946.820
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge, data = repro_nfC, family = "gaussian", na.action = na.fail)
model_selection = dredge(FC1)
model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge,
## data = repro_nfC, family = "gaussian", na.action = na.fail,
## ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 2 -22.74 + 1.305 3 -966.679
## 4 -24.62 + 1.311 + 7 -963.982
## 6 -22.39 + 1.295 + 5 -966.198
## 8 -24.57 + 1.305 + + 9 -963.461
## 16 -24.28 + 1.314 + + + 15 -962.257
## 5 47.93 + + 4 -1078.562
## 1 48.49 + 2 -1083.107
## 7 50.74 + + + 8 -1077.345
## 3 52.70 + + 6 -1079.491
## 15 46.80 + + + + 14 -1076.437
## AICc delta weight
## 2 1939.5 0.00 0.682
## 4 1942.5 3.01 0.151
## 6 1942.7 3.20 0.137
## 8 1945.7 6.28 0.029
## 16 1956.8 17.31 0.000
## 5 2165.3 225.84 0.000
## 1 2170.3 230.80 0.000
## 7 2171.3 231.88 0.000
## 3 2171.4 231.90 0.000
## 15 2182.8 243.38 0.000
## Models ranked by AICc(x)
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.776 0.628 0.604 0.564 0.528 0.548 0.648 0.632 0.548 0.436 0.528 0.432 0.496 0.436 0.432 0.452 0.724 0.516 0.784 0.664 ...
selected_models = get.models(model_selection, 1:nrow(model_selection))
tab_model(selected_models, show.intercept = F) # No effect
| Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | |||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| DeathAge | 1.30 | 1.18 – 1.43 | <0.001 | 1.31 | 1.18 – 1.44 | <0.001 | 1.29 | 1.17 – 1.42 | <0.001 | 1.31 | 1.18 – 1.43 | <0.001 | 1.31 | 1.18 – 1.44 | <0.001 | |||||||||||||||
| Isolate [C14] | 4.03 | -2.97 – 11.04 | 0.259 | 3.83 | -3.22 – 10.88 | 0.287 | 3.21 | -4.95 – 11.37 | 0.441 | 1.19 | -10.42 – 12.80 | 0.841 | 0.82 | -10.79 – 12.43 | 0.890 | 5.39 | -8.07 – 18.85 | 0.433 | ||||||||||||
| Isolate [C18] | 3.72 | -3.09 – 10.52 | 0.284 | 3.82 | -3.00 – 10.63 | 0.272 | 3.17 | -4.64 – 10.98 | 0.426 | -3.36 | -14.52 – 7.80 | 0.555 | -4.35 | -15.56 – 6.87 | 0.447 | 0.98 | -11.90 – 13.86 | 0.882 | ||||||||||||
| Isolate [C24] | -2.42 | -9.45 – 4.62 | 0.501 | -2.51 | -9.64 – 4.63 | 0.491 | -2.96 | -11.08 – 5.16 | 0.475 | -5.37 | -17.12 – 6.37 | 0.370 | -6.47 | -18.12 – 5.19 | 0.277 | -1.23 | -14.62 – 12.17 | 0.858 | ||||||||||||
| Isolate [CTRL] | 1.21 | -7.28 – 9.69 | 0.780 | 2.25 | -7.51 – 12.00 | 0.652 | -0.81 | -18.17 – 16.54 | 0.927 | -7.95 | -23.93 – 8.03 | 0.330 | -15.55 | -29.36 – -1.74 | 0.027 | 7.65 | -20.95 – 36.26 | 0.600 | ||||||||||||
| Temperature [linear] | 1.63 | -2.10 – 5.37 | 0.391 | 1.71 | -2.54 – 5.96 | 0.431 | -1.60 | -14.45 – 11.26 | 0.808 | 9.22 | 3.23 – 15.20 | 0.003 | 7.35 | 0.41 – 14.29 | 0.038 | 19.74 | -1.19 – 40.67 | 0.065 | ||||||||||||
| Temperature [quadratic] | 0.81 | -3.19 – 4.80 | 0.692 | 0.98 | -3.23 – 5.18 | 0.649 | 1.01 | -8.88 – 10.90 | 0.841 | -2.19 | -8.71 – 4.33 | 0.510 | -1.09 | -8.01 – 5.83 | 0.757 | -8.20 | -24.46 – 8.05 | 0.323 | ||||||||||||
|
Isolate [C14] × Temperature [linear] |
3.67 | -11.49 – 18.83 | 0.635 | -14.13 | -38.99 – 10.72 | 0.265 | ||||||||||||||||||||||||
|
Isolate [C18] × Temperature [linear] |
3.79 | -11.12 – 18.70 | 0.618 | -14.78 | -39.20 – 9.63 | 0.235 | ||||||||||||||||||||||||
|
Isolate [C24] × Temperature [linear] |
3.84 | -10.96 – 18.64 | 0.611 | -12.61 | -36.88 – 11.67 | 0.309 | ||||||||||||||||||||||||
|
Isolate [C14] × Temperature [quadratic] |
-2.63 | -15.81 – 10.55 | 0.696 | 8.53 | -13.14 – 30.20 | 0.440 | ||||||||||||||||||||||||
|
Isolate [C18] × Temperature [quadratic] |
4.00 | -8.12 – 16.11 | 0.518 | 9.27 | -10.71 – 29.24 | 0.363 | ||||||||||||||||||||||||
|
Isolate [C24] × Temperature [quadratic] |
-2.94 | -16.34 – 10.46 | 0.667 | 5.90 | -16.16 – 27.96 | 0.600 | ||||||||||||||||||||||||
| Observations | 228 | 228 | 228 | 228 | 228 | 228 | 228 | 228 | 228 | 228 | ||||||||||||||||||||
| R2 / R2 adjusted | 0.640 / 0.637 | 0.648 / 0.639 | 0.641 / 0.635 | 0.650 / 0.637 | 0.654 / 0.631 | 0.039 / 0.026 | 0.000 / -0.004 | 0.049 / 0.019 | 0.031 / 0.009 | 0.057 / -0.000 | ||||||||||||||||||||
summary(best_model)
## Family: gaussian ( identity )
## Formula: Nb_offspring ~ DeathAge + 1
## Data: repro_nfC
##
## AIC BIC logLik deviance df.resid
## 1939.4 1949.6 -966.7 1933.4 225
##
##
## Dispersion estimate for gaussian family (sigma^2): 282
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.74119 3.70948 -6.131 8.76e-10 ***
## DeathAge 1.30473 0.06482 20.127 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Again, the husbandry control do not change the conclusions.
Here we just described the average number of offspring produced per day in total in each group
# Function to compute mean and 95% confidence interval
mean_ci = function(x) {
x = x[!is.na(x)]
m = mean(x)
se = sd(x) / sqrt(length(x))
c(mean = m,
CI_low = m - 1.96 * se,
CI_high = m + 1.96 * se)
}
# Define groups
g1 = mean_ci(repro_nfC$repro_rate[repro_nfC$Isolate == "CTRL"])
g2 = mean_ci(repro_nf$repro_rate[repro_nfC$Temperature == 20])
g3 = mean_ci(repro_nf$repro_rate)
g4 = mean_ci(repro_inf$repro_rate)
rbind(
noninfected_unexposed = g1,
noninfected_exposed_20 = g2,
noninfected_exposed_all = g3,
infected_all = g4
)
## mean CI_low CI_high
## noninfected_unexposed 0.7203259 0.5747725 0.8658793
## noninfected_exposed_20 0.7463513 0.6490932 0.8436094
## noninfected_exposed_all 0.8219201 0.7675482 0.8762921
## infected_all 0.1699311 0.1302886 0.2095736
This section explores how temperature and the genotype of the parasite contribute to the abundance of each stage of the life-cycle after infection, using generalized linear mixed models.
Cauliflowers corresponds to the first stage of the parasite life-cycle and here we are interested to understand how their presence/absence is affected by the temperature and/or the genotype of the spores. We also investigated the potential impact of the variations in time post-exposure, as all screened individuals died at different time.
data_c1 = subset(data, data$cauliflowers != "not-well-counted" & data$cauliflowers != "counted-after-freezer" & data$stade_max != "2nd_Cauliflower" & Temperature != "NA")
data_c1$Temperature = factor(data_c1$Temperature, ordered = TRUE)
Caul = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate+DeathAge, family = "binomial", REML= T)
CaulR = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate+DeathAge + (1|Batch), family = "binomial", REML= T)
AICc(Caul, CaulR)
## df AICc
## Caul 7 119.6534
## CaulR 8 121.7682
Caul1 = glmmTMB(data=data_c1, cauliflowers~Temperature*Isolate*DeathAge, family = "binomial", REML= F)
Caul2 = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate*DeathAge, family = "binomial", REML= F)
Caul3 = glmmTMB(data=data_c1, cauliflowers~Temperature*Isolate+DeathAge, family = "binomial", REML= F)
Caul4 = glmmTMB(data=data_c1, cauliflowers~Isolate+DeathAge*Temperature, family = "binomial", REML= F)
Caul5 = glmmTMB(data=data_c1, cauliflowers~Isolate*DeathAge, family = "binomial", REML= F)
Caul6 = glmmTMB(data=data_c1, cauliflowers~Temperature*Isolate, family = "binomial", REML= F)
Caul7 = glmmTMB(data=data_c1, cauliflowers~DeathAge*Temperature, family = "binomial", REML= F)
Caul8 = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate+DeathAge, family = "binomial", REML= F)
Caul9 = glmmTMB(data=data_c1, cauliflowers~Isolate+DeathAge, family = "binomial", REML= F)
Caul10 = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate, family = "binomial", REML= F)
Caul11 = glmmTMB(data=data_c1, cauliflowers~DeathAge+Temperature, family = "binomial", REML= F)
Caul12 = glmmTMB(data=data_c1, cauliflowers~DeathAge, family = "binomial", REML= F)
Caul13 = glmmTMB(data=data_c1, cauliflowers~Isolate, family = "binomial", REML= F)
Caul14 = glmmTMB(data=data_c1, cauliflowers~Temperature, family = "binomial", REML= F)
Caul0 = glmmTMB(data=data_c1, cauliflowers~1, family = "binomial", REML= F)
AICc(Caul1, Caul2, Caul3, Caul4, Caul5, Caul6, Caul7, Caul8, Caul9, Caul10, Caul11, Caul12, Caul13, Caul14, Caul0)
## df AICc
## Caul1 24 NA
## Caul2 10 121.1953
## Caul3 13 123.4742
## Caul4 9 121.1437
## Caul5 8 121.2525
## Caul6 12 151.2119
## Caul7 6 115.2956
## Caul8 7 116.9950
## Caul9 5 117.3381
## Caul10 6 148.7161
## Caul11 4 111.2346
## Caul12 2 111.9709
## Caul13 4 157.5159
## Caul14 3 144.1507
## Caul0 1 152.8702
tab_model(Caul11, Caul12, Caul7, Caul8, Caul5, Caul2, Caul4, show.intercept = F)
| cauliflowers | cauliflowers | cauliflowers | cauliflowers | cauliflowers | cauliflowers | cauliflowers | |||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| DeathAge | 0.91 | 0.87 – 0.95 | <0.001 | 0.90 | 0.87 – 0.94 | <0.001 | 0.91 | 0.87 – 0.95 | <0.001 | 0.91 | 0.87 – 0.95 | <0.001 | 0.85 | 0.73 – 1.00 | 0.049 | 0.86 | 0.73 – 1.01 | 0.064 | 0.91 | 0.87 – 0.95 | <0.001 |
| Temperature [linear] | 0.34 | 0.11 – 1.05 | 0.060 | 0.24 | 0.02 – 3.09 | 0.271 | 0.34 | 0.11 – 1.08 | 0.066 | 0.34 | 0.10 – 1.08 | 0.068 | 0.23 | 0.02 – 3.33 | 0.284 | ||||||
| Temperature [quadratic] | 0.97 | 0.36 – 2.57 | 0.944 | 0.84 | 0.08 – 9.24 | 0.886 | 0.96 | 0.36 – 2.57 | 0.940 | 0.94 | 0.35 – 2.51 | 0.895 | 0.81 | 0.07 – 9.22 | 0.866 | ||||||
|
DeathAge × Temperature [linear] |
1.01 | 0.93 – 1.10 | 0.753 | 1.01 | 0.93 – 1.10 | 0.755 | |||||||||||||||
|
DeathAge × Temperature [quadratic] |
1.00 | 0.92 – 1.09 | 0.911 | 1.01 | 0.92 – 1.10 | 0.890 | |||||||||||||||
| Isolate [C1] | 0.91 | 0.21 – 3.99 | 0.897 | 0.16 | 0.00 – 12.00 | 0.409 | 0.11 | 0.00 – 9.38 | 0.331 | 0.89 | 0.20 – 3.94 | 0.878 | |||||||||
| Isolate [C18] | 0.85 | 0.19 – 3.83 | 0.831 | 0.06 | 0.00 – 4.98 | 0.214 | 0.10 | 0.00 – 10.12 | 0.331 | 0.87 | 0.19 – 3.99 | 0.857 | |||||||||
| Isolate [C24] | 1.33 | 0.34 – 5.16 | 0.683 | 0.91 | 0.01 – 134.04 | 0.969 | 1.02 | 0.01 – 159.66 | 0.994 | 1.34 | 0.34 – 5.23 | 0.676 | |||||||||
| Isolate [C1] × DeathAge | 1.08 | 0.91 – 1.29 | 0.389 | 1.09 | 0.91 – 1.31 | 0.340 | |||||||||||||||
| Isolate [C18] × DeathAge | 1.10 | 0.92 – 1.31 | 0.285 | 1.09 | 0.91 – 1.30 | 0.363 | |||||||||||||||
| Isolate [C24] × DeathAge | 1.02 | 0.82 – 1.25 | 0.886 | 1.01 | 0.82 – 1.25 | 0.899 | |||||||||||||||
| Observations | 291 | 291 | 291 | 291 | 291 | 291 | 291 | ||||||||||||||
| R2 Tjur | 0.194 | 0.154 | 0.197 | 0.203 | 0.183 | 0.220 | 0.205 | ||||||||||||||
# No differences
simulateResiduals(Caul11, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3388901 0.3332647 0.415397 0.04419869 0.3642384 0.5669417 0.7992067 0.6843859 0.5363047 0.5255567 0.4620051 0.9983571 0.2413181 0.1928017 0.1345064 0.2465995 0.7431666 0.4926976 0.5356813 0.0344749 ...
# best model
Caul11 = glmmTMB(data=data_c1, cauliflowers~DeathAge+Temperature, family = "binomial", REML=T)
summary(Caul11)
## Family: binomial ( logit )
## Formula: cauliflowers ~ DeathAge + Temperature
## Data: data_c1
##
## AIC BIC logLik deviance df.resid
## 116.8 131.4 -54.4 108.8 287
##
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.66549 0.67120 0.991 0.3214
## DeathAge -0.09664 0.02279 -4.241 2.23e-05 ***
## Temperature.L -1.08254 0.57628 -1.878 0.0603 .
## Temperature.Q -0.03522 0.49967 -0.070 0.9438
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(Caul11, show.intercept = F)
| cauliflowers | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| DeathAge | 0.91 | 0.87 – 0.95 | <0.001 |
| Temperature [linear] | 0.34 | 0.11 – 1.05 | 0.060 |
| Temperature [quadratic] | 0.97 | 0.36 – 2.57 | 0.944 |
| Observations | 291 | ||
| R2 Tjur | 0.194 | ||
We observed that the time post-exposure is reducing of 10% the
probability to observe first stage cauliflowers, temperature might have
a role that is not significant with the current effective of our
experimentation.
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_line()`).
The grape seed stage corresponds to the second visible stage of the parasite life cycle and here we are interested in understanding how its presence/absence is affected by temperature and the genotype of the spores. We also investigated the potential impact of variations in time post-exposure, as all screened individuals died at different times.
data_crushed_inf_c = subset(infected, infected$grapes!="counted-after-freezer" & infected$grapes!="counted-after-freezer")
Grap = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate+DeathAge, family = "binomial", REML=F)
GrapR = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate+DeathAge+ (1|Batch), family = "binomial", REML=F)
AICc(Grap,GrapR) # No need for random effect
## df AICc
## Grap 7 209.1698
## GrapR 8 210.7109
Grap1 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature*Isolate*DeathAge, family = "binomial")
Grap2 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate*DeathAge, family = "binomial")
Grap3 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature*Isolate+DeathAge, family = "binomial")
Grap4 = glmmTMB(data=data_crushed_inf_c, grapes~Isolate+DeathAge*Temperature, family = "binomial")
Grap5 = glmmTMB(data=data_crushed_inf_c, grapes~Isolate*DeathAge, family = "binomial")
Grap6 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature*Isolate, family = "binomial")
Grap7 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature*DeathAge, family = "binomial")
Grap8 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate+DeathAge, family = "binomial")
Grap9 = glmmTMB(data=data_crushed_inf_c, grapes~Isolate+DeathAge, family = "binomial")
Grap10 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate, family = "binomial")
Grap11 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+DeathAge, family = "binomial")
Grap12 = glmmTMB(data=data_crushed_inf_c, grapes~DeathAge, family = "binomial")
Grap13 = glmmTMB(data=data_crushed_inf_c, grapes~Isolate, family = "binomial")
Grap14 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature, family = "binomial")
Grap0 = glmmTMB(data=data_crushed_inf_c, grapes~1, family = "binomial")
AICc(Grap1, Grap2, Grap3, Grap4, Grap5, Grap6, Grap7, Grap8, Grap9, Grap10, Grap11, Grap12, Grap13, Grap14, Grap0)
## df AICc
## Grap1 24 NA
## Grap2 10 204.7698
## Grap3 13 219.5599
## Grap4 9 212.0557
## Grap5 8 207.6257
## Grap6 12 217.6432
## Grap7 6 206.7304
## Grap8 7 209.1698
## Grap9 5 213.9978
## Grap10 6 207.8652
## Grap11 4 204.1768
## Grap12 2 209.5133
## Grap13 4 213.0738
## Grap14 3 202.8924
## Grap0 1 208.6435
tab_model(Grap14, Grap11, Grap10, Grap8, Grap0, show.intercept = F)
| grapes | grapes | grapes | grapes | grapes | |||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| Temperature [linear] | 2.69 | 1.38 – 5.25 | 0.004 | 2.66 | 1.36 – 5.20 | 0.004 | 2.70 | 1.37 – 5.32 | 0.004 | 2.67 | 1.35 – 5.27 | 0.005 | |||
| Temperature [quadratic] | 1.44 | 0.79 – 2.62 | 0.237 | 1.41 | 0.77 – 2.57 | 0.268 | 1.40 | 0.76 – 2.56 | 0.278 | 1.37 | 0.74 – 2.51 | 0.315 | |||
| DeathAge | 1.01 | 0.99 – 1.04 | 0.366 | 1.01 | 0.99 – 1.04 | 0.348 | |||||||||
| Isolate [C1] | 1.76 | 0.63 – 4.89 | 0.281 | 1.84 | 0.65 – 5.18 | 0.249 | |||||||||
| Isolate [C18] | 1.57 | 0.59 – 4.17 | 0.363 | 1.44 | 0.53 – 3.90 | 0.478 | |||||||||
| Isolate [C24] | 1.25 | 0.47 – 3.36 | 0.652 | 1.22 | 0.45 – 3.29 | 0.693 | |||||||||
| Observations | 158 | 158 | 158 | 158 | 158 | ||||||||||
| R2 Tjur | 0.058 | 0.064 | 0.068 | 0.075 | 0.000 | ||||||||||
# Temperature increase the chances to observe grapes
simulateResiduals(Grap14, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.7676696 0.05811451 0.2656332 0.3575351 0.9665073 0.3685181 0.7259261 0.04275554 0.06533833 0.1598172 0.02149855 0.8880048 0.4096685 0.3173512 0.08902606 0.6813057 0.1156503 0.1831655 0.3961829 0.6764339 ...
summary(Grap14)
## Family: binomial ( logit )
## Formula: grapes ~ Temperature
## Data: data_crushed_inf_c
##
## AIC BIC logLik deviance df.resid
## 202.7 211.9 -98.4 196.7 155
##
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7337 0.1873 3.918 8.94e-05 ***
## Temperature.L 0.9886 0.3416 2.894 0.0038 **
## Temperature.Q 0.3618 0.3062 1.182 0.2373
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(Grap14, show.intercept = F)
| grapes | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| Temperature [linear] | 2.69 | 1.38 – 5.25 | 0.004 |
| Temperature [quadratic] | 1.44 | 0.79 – 2.62 | 0.237 |
| Observations | 158 | ||
| R2 Tjur | 0.058 | ||
We observed that the temperature increases the odds of observing grape
stage seeds by 168%.
After grape seed formation, spores begin their final maturation stage and develop their final shape. However, this is gradual and when an individual died we could observe partially matured spores, which we called immature spores. Here we are interested in understanding how their presence/absence is affected by temperature and/or the genotype of the spores. We also investigated the potential impact of variations in time post-exposure, as all screened individuals died at different times.
data_crushed_inf_c = subset(infected, infected$imm_load!="counted-after-freezer")
data_crushed_inf_c$imm_presence = ifelse(data_crushed_inf_c$imm >0, 1,0)
ImmP = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate+DeathAge, family = "binomial")
ImmPR = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate+DeathAge+ (1|Batch), family = "binomial")
AICc(ImmP,ImmPR) # No need for random effect
## df AICc
## ImmP 7 109.2572
## ImmPR 8 108.0838
ImmP1 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature*Isolate*DeathAge, family = "binomial")
ImmP2 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate*DeathAge, family = "binomial")
ImmP3 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature*Isolate+DeathAge, family = "binomial")
ImmP4 = glmmTMB(data=data_crushed_inf_c, imm_presence~Isolate+DeathAge*Temperature, family = "binomial")
ImmP5 = glmmTMB(data=data_crushed_inf_c, imm_presence~Isolate*DeathAge, family = "binomial")
ImmP6 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature*Isolate, family = "binomial")
ImmP7 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature*DeathAge, family = "binomial")
ImmP8 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate+DeathAge, family = "binomial")
ImmP9 = glmmTMB(data=data_crushed_inf_c, imm_presence~Isolate+DeathAge, family = "binomial")
ImmP10 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate, family = "binomial")
ImmP11 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+DeathAge, family = "binomial")
ImmP12 = glmmTMB(data=data_crushed_inf_c, imm_presence~DeathAge, family = "binomial")
ImmP13 = glmmTMB(data=data_crushed_inf_c, imm_presence~Isolate, family = "binomial")
ImmP14 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature, family = "binomial")
ImmP0 = glmmTMB(data=data_crushed_inf_c, imm_presence~1, family = "binomial")
AICc(ImmP1, ImmP2, ImmP3, ImmP4, ImmP5, ImmP6, ImmP7, ImmP8, ImmP9, ImmP10, ImmP11, ImmP12, ImmP13, ImmP14, ImmP0)
## df AICc
## ImmP1 24 NA
## ImmP2 10 111.6544
## ImmP3 13 112.0329
## ImmP4 9 110.5120
## ImmP5 8 110.1232
## ImmP6 12 146.8456
## ImmP7 6 104.4657
## ImmP8 7 109.2572
## ImmP9 5 108.3799
## ImmP10 6 148.4499
## ImmP11 4 103.2475
## ImmP12 2 102.5769
## ImmP13 4 147.5806
## ImmP14 3 143.6712
## ImmP0 1 143.0355
tab_model(ImmP12, ImmP11, ImmP7, show.intercept = F)
| imm_presence | imm_presence | imm_presence | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| DeathAge | 1.13 | 1.08 – 1.19 | <0.001 | 1.14 | 1.08 – 1.19 | <0.001 | 1.13 | 1.07 – 1.20 | <0.001 |
| Temperature [linear] | 1.20 | 0.48 – 3.00 | 0.690 | 1.48 | 0.03 – 80.40 | 0.849 | |||
| Temperature [quadratic] | 0.43 | 0.15 – 1.20 | 0.106 | 0.04 | 0.00 – 0.99 | 0.050 | |||
|
Temperature [linear] × DeathAge |
0.99 | 0.89 – 1.10 | 0.882 | ||||||
|
Temperature [quadratic] × DeathAge |
1.08 | 0.99 – 1.17 | 0.101 | ||||||
| Observations | 167 | 167 | 167 | ||||||
| R2 Tjur | 0.311 | 0.344 | 0.361 | ||||||
# Temperature increase the chances to observe imm_presence
simulateResiduals(ImmP12, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.08681092 0.4781304 0.5053635 0.3460266 0.9493186 0.5617266 0.07295504 0.2366983 0.1756649 0.6546839 0.7908228 0.8876897 0.893821 0.7406902 0.239505 0.7249811 0.1957589 0.5392972 0.5591017 0.8727609 ...
summary(ImmP12)
## Family: binomial ( logit )
## Formula: imm_presence ~ DeathAge
## Data: data_crushed_inf_c
##
## AIC BIC logLik deviance df.resid
## 102.5 108.7 -49.3 98.5 165
##
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.12919 0.84637 -3.697 0.000218 ***
## DeathAge 0.12522 0.02326 5.384 7.27e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(ImmP12, show.intercept = F)
| imm_presence | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| DeathAge | 1.13 | 1.08 – 1.19 | <0.001 |
| Observations | 167 | ||
| R2 Tjur | 0.311 | ||
We observed that only the time between exposure and death is increase of
13% the probability to observe immmature presence after host death,
which is why we removed this abalysis from the main manuscript.
Spores reach their final maturation stage and can now infect new hosts. We are interested in understanding how individual mature parasite load is affected by temperature and/or the genotype of the spores. We also investigated the potential impact of variations in time post-exposure, as all screened individuals died at different times.
data_mat = subset(infected, infected$mature_load>0)
data_mat$Temperature = factor(data_mat$Temperature, ordered = F)
mat1 = glmmTMB(data=data_mat, mature_load/10000000~Temperature*Isolate*DeathAge, family = "gaussian", REML=T)
mat1R = glmmTMB(data=data_mat, mature_load/10000000~Temperature*Isolate*DeathAge+(1|Batch), family = "gaussian", REML=T)
AICc(mat1, mat1R) # No need of a batch effect
## df AICc
## mat1 25 242.8509
## mat1R 26 245.6588
mat2 = glmmTMB(data=data_mat, mature_load/10000000~Temperature+Isolate*DeathAge, family = "gaussian")
mat3 = glmmTMB(data=data_mat, mature_load/10000000~Temperature*Isolate+DeathAge, family = "gaussian")
mat4 = glmmTMB(data=data_mat, mature_load/10000000~Isolate+DeathAge*Temperature, family = "gaussian")
mat5 = glmmTMB(data=data_mat, mature_load/10000000~Isolate*DeathAge, family = "gaussian")
mat6 = glmmTMB(data=data_mat, mature_load/10000000~Temperature*Isolate, family = "gaussian")
mat7 = glmmTMB(data=data_mat, mature_load/10000000~Temperature*DeathAge, family = "gaussian")
mat8 = glmmTMB(data=data_mat, mature_load/10000000~Temperature+Isolate+DeathAge, family = "gaussian")
mat9 = glmmTMB(data=data_mat, mature_load/10000000~Isolate+DeathAge, family = "gaussian")
mat10 = glmmTMB(data=data_mat, mature_load/10000000~Temperature+Isolate, family = "gaussian")
mat11 = glmmTMB(data=data_mat, mature_load/10000000~Temperature+DeathAge, family = "gaussian")
mat12 = glmmTMB(data=data_mat, mature_load/10000000~DeathAge, family = "gaussian")
mat13 = glmmTMB(data=data_mat, mature_load/10000000~Isolate, family = "gaussian")
mat14 = glmmTMB(data=data_mat, mature_load/10000000~Temperature, family = "gaussian")
mat0 = glmmTMB(data=data_mat, mature_load/10000000~1, family = "gaussian")
AICc(mat1, mat2, mat3, mat4, mat5, mat6, mat7, mat8, mat9, mat10, mat11, mat12, mat13, mat14, mat0)
## df AICc
## mat1 25 242.85092
## mat2 11 97.42244
## mat3 14 104.23161
## mat4 10 99.48683
## mat5 9 93.04474
## mat6 13 150.21483
## mat7 7 107.51440
## mat8 8 96.60519
## mat9 6 92.28979
## mat10 7 139.42053
## mat11 5 104.35369
## mat12 3 100.24447
## mat13 5 136.94938
## mat14 4 154.88560
## mat0 2 152.62650
tab_model(mat9, mat5, show.intercept = F)
| mature_load/1e+07 | mature_load/1e+07 | |||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p |
| Isolate [C1] | 0.09 | -0.07 – 0.24 | 0.286 | -0.70 | -1.46 – 0.06 | 0.071 |
| Isolate [C18] | 0.28 | 0.13 – 0.43 | <0.001 | -0.44 | -1.10 – 0.22 | 0.189 |
| Isolate [C24] | 0.09 | -0.07 – 0.24 | 0.286 | -0.35 | -1.03 – 0.33 | 0.308 |
| DeathAge | 0.02 | 0.01 – 0.02 | <0.001 | 0.01 | -0.00 – 0.02 | 0.122 |
| Isolate [C1] × DeathAge | 0.02 | 0.00 – 0.03 | 0.043 | |||
| Isolate [C18] × DeathAge | 0.01 | 0.00 – 0.03 | 0.028 | |||
| Isolate [C24] × DeathAge | 0.01 | -0.00 – 0.02 | 0.184 | |||
| Observations | 145 | 145 | ||||
| R2 / R2 adjusted | 0.378 / 0.356 | 0.403 / 0.368 | ||||
simulateResiduals(mat5, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.468 0.368 0.34 0.828 0.316 0.808 0.128 0.4 0.48 0.48 0.124 0.096 0.224 0.332 0.36 0.516 0.54 0.3 0.34 0.48 ...
tab_model(mat5, show.intercept = F)
| mature_load/1e+07 | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Isolate [C1] | -0.70 | -1.46 – 0.06 | 0.071 |
| Isolate [C18] | -0.44 | -1.10 – 0.22 | 0.189 |
| Isolate [C24] | -0.35 | -1.03 – 0.33 | 0.308 |
| DeathAge | 0.01 | -0.00 – 0.02 | 0.122 |
| Isolate [C1] × DeathAge | 0.02 | 0.00 – 0.03 | 0.043 |
| Isolate [C18] × DeathAge | 0.01 | 0.00 – 0.03 | 0.028 |
| Isolate [C24] × DeathAge | 0.01 | -0.00 – 0.02 | 0.184 |
| Observations | 145 | ||
| R2 / R2 adjusted | 0.403 / 0.368 | ||
C1 and C18 have significantly increased growth speed over time compared
to C14 (référence) and C24.
Cauliflowers can also be observed if parasites restart a second cycle inside the same host. Here we are interested in understanding how their presence/absence is affected by temperature and the genotype of the spores. We also investigated the potential impact of variations in time post-exposure, as all screened individuals died at different times.
data_c2 = subset(infected, infected$cauliflowers!="counted-after-freezer" & infected$cauliflowers!="counted-after-freezer"& infected$stade_max!="Cauliflower")
data_c2$Temperature = factor(data_c2$Temperature, ordered = T)
Caul2 = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate+DeathAge, family = "binomial")
Caul2R = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate+DeathAge + (1|Batch), family = "binomial", REML= T)
AICc(Caul2, Caul2R)
## df AICc
## Caul2 13 207.2512
## Caul2R 14 207.1492
Caul2.1 = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate*DeathAge, family = "binomial")
Caul2.2 = glmmTMB(data=data_c2, cauliflowers~Temperature+Isolate*DeathAge, family = "binomial")
Caul2.3 = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate+DeathAge, family = "binomial")
Caul2.4 = glmmTMB(data=data_c2, cauliflowers~Isolate+DeathAge*Temperature, family = "binomial")
Caul2.5 = glmmTMB(data=data_c2, cauliflowers~Isolate*DeathAge, family = "binomial")
Caul2.6 = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate, family = "binomial")
Caul2.7 = glmmTMB(data=data_c2, cauliflowers~DeathAge*Temperature, family = "binomial")
Caul2.8 = glmmTMB(data=data_c2, cauliflowers~Temperature+Isolate+DeathAge, family = "binomial")
Caul2.9 = glmmTMB(data=data_c2, cauliflowers~Isolate+DeathAge, family = "binomial")
Caul2.10 = glmmTMB(data=data_c2, cauliflowers~Temperature+Isolate, family = "binomial")
Caul2.11 = glmmTMB(data=data_c2, cauliflowers~DeathAge+Temperature, family = "binomial")
Caul2.12 = glmmTMB(data=data_c2, cauliflowers~DeathAge, family = "binomial")
Caul2.13 = glmmTMB(data=data_c2, cauliflowers~Isolate, family = "binomial")
Caul2.14 = glmmTMB(data=data_c2, cauliflowers~Temperature, family = "binomial")
Caul2.0 = glmmTMB(data=data_c2, cauliflowers~1, family = "binomial")
AICc(Caul2.1, Caul2.2, Caul2.3, Caul2.4, Caul2.5, Caul2.6, Caul2.7, Caul2.8, Caul2.9, Caul2.10, Caul2.11, Caul2.12, Caul2.13, Caul2.14, Caul2.0)
## df AICc
## Caul2.1 24 221.5962
## Caul2.2 10 206.3934
## Caul2.3 13 207.2512
## Caul2.4 9 206.1231
## Caul2.5 8 205.8774
## Caul2.6 12 205.5078
## Caul2.7 6 205.5147
## Caul2.8 7 201.9549
## Caul2.9 5 201.6345
## Caul2.10 6 201.1274
## Caul2.11 4 201.8077
## Caul2.12 2 201.3683
## Caul2.13 4 200.2414
## Caul2.14 3 199.9651
## Caul2.0 1 199.3769
tab_model(Caul2.0, Caul2.14, Caul2.13, Caul2.12, Caul2.10, show.intercept = F)
| cauliflowers | cauliflowers | cauliflowers | cauliflowers | cauliflowers | |||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| Temperature [linear] | 1.06 | 0.60 – 1.89 | 0.837 | 1.04 | 0.58 – 1.89 | 0.891 | |||||||||
| Temperature [quadratic] | 0.58 | 0.32 – 1.05 | 0.070 | 0.57 | 0.31 – 1.05 | 0.073 | |||||||||
| Isolate [C1] | 0.91 | 0.32 – 2.61 | 0.858 | 0.98 | 0.33 – 2.85 | 0.965 | |||||||||
| Isolate [C18] | 0.38 | 0.14 – 1.02 | 0.054 | 0.39 | 0.14 – 1.07 | 0.067 | |||||||||
| Isolate [C24] | 0.69 | 0.25 – 1.94 | 0.488 | 0.72 | 0.25 – 2.05 | 0.535 | |||||||||
| DeathAge | 1.00 | 0.98 – 1.03 | 0.799 | ||||||||||||
| Observations | 145 | 145 | 145 | 145 | 145 | ||||||||||
| R2 Tjur | 0.000 | 0.024 | 0.037 | 0.000 | 0.059 | ||||||||||
# No differences
simulateResiduals(Caul2.13, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2323121 0.4716365 0.7126411 0.6861518 0.2546818 0.3056335 0.645681 0.6405459 0.9225231 0.8999696 0.9554663 0.1963001 0.20196 0.6420615 0.6905008 0.00794358 0.1867782 0.665696 0.4064166 0.2371338 ...
tab_model(Caul2.10, show.intercept = F)
| cauliflowers | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| Temperature [linear] | 1.04 | 0.58 – 1.89 | 0.891 |
| Temperature [quadratic] | 0.57 | 0.31 – 1.05 | 0.073 |
| Isolate [C1] | 0.98 | 0.33 – 2.85 | 0.965 |
| Isolate [C18] | 0.39 | 0.14 – 1.07 | 0.067 |
| Isolate [C24] | 0.72 | 0.25 – 2.05 | 0.535 |
| Observations | 145 | ||
| R2 Tjur | 0.059 | ||
We observed that the null model is among the best ones and that the
measured effect are relatively low, so we preffered the model taking
into account the effect of the genotype on the second cauliflower
formation, which seem more biologically relevant, however it should be
interptreted with caution.
The isolate C18 is a little bit less likely to start a second cycle
inside it’s host thant the other lineages (66% less in average).
The purpose of this section is to analyze the impact of temperature and parasite genotype on the maximum stage reached by the parasite at host death. We also included the impact of host age at death as individuals died at different times post-exposure.
infected = subset(infected,stade_max!="None")
infected$stade_max = factor(infected$stade_max, ordered = T)
infected$Temperature = factor(infected$Temperature, ordered = T)
# On crée une nouvelle variable pour ne pas écraser l'originale
infected$stade_simplifie <- fct_collapse(infected$stade_max,
"Caul_Grapes" = c("Cauliflower", "Grapes"),
"Imm_Mature" = c("Immature", "Mature"),
"2nd_Caul" = "2nd_Cauliflower"
)
# Crucial : On s'assure que R comprend que c'est ordonné
# (Vérifie bien que cet ordre est le bon ordre biologique !)
infected$stade_simplifie <- factor(infected$stade_simplifie,
levels = c("Caul_Grapes", "Imm_Mature", "2nd_Caul"),
ordered = TRUE)
#modele13 = polr(stade_max ~ Temperature * Isolate * DeathAge, data = infected, Hess = TRUE)
modele14 = polr(stade_max ~ Temperature + Isolate + DeathAge, data = infected, Hess = TRUE)
modele14R = clmm(stade_max ~ Temperature + Isolate + DeathAge + (1 | Batch), data = infected)
AICc(modele14R, modele14) # No need for random factor
## df AICc
## modele14R 11 293.3242
## modele14 10 292.0298
#modele1 = polr(stade_max ~ Temperature + Isolate * DeathAge, data = infected, Hess = TRUE)
modele2 = polr(stade_max ~ Temperature * Isolate + DeathAge, data = infected, Hess = TRUE)
modele3 = polr(stade_max ~ Isolate + DeathAge * Temperature, data = infected, Hess = TRUE)
#modele4 = polr(stade_max ~ Isolate * DeathAge, data = infected, Hess = TRUE)
modele5 = polr(stade_max ~ Temperature * Isolate, data = infected, Hess = TRUE)
modele6 = polr(stade_max ~ DeathAge * Temperature, data = infected, Hess = TRUE)
modele7 = polr(stade_max ~ Isolate + DeathAge, data = infected, Hess = TRUE)
modele8 = polr(stade_max ~ Temperature + Isolate, data = infected, Hess = TRUE)
modele9 = polr(stade_max ~ DeathAge + Temperature, data = infected, Hess = TRUE)
modele10 = polr(stade_max ~ Temperature, data = infected, Hess = TRUE)
modele11 = polr(stade_max ~ Isolate, data = infected, Hess = TRUE)
modele12 = polr(stade_max ~ DeathAge, data = infected, Hess = TRUE)
modele0 = polr(stade_max ~ 1, data = infected, Hess = TRUE)
AICc(modele2,modele3,modele5,modele6,modele7,modele8,modele9,modele10,modele11,modele12,modele14,modele0)
## df AICc
## modele2 16 301.0434
## modele3 12 293.9404
## modele5 15 347.5868
## modele6 9 297.1830
## modele7 8 299.5568
## modele8 9 341.0864
## modele9 7 294.9343
## modele10 6 335.5819
## modele11 7 342.4234
## modele12 5 301.1638
## modele14 10 292.0298
## modele0 4 337.1707
tab_model(modele14, modele3, modele9)
| stade_max | stade_max | stade_max | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| Cauliflower|Grapes | 3.24 | 0.85 – 12.35 | 0.084 | 2.01 | 0.44 – 9.15 | 0.364 | 2.60 | 0.81 – 8.37 | 0.109 |
| Grapes|Immature | 3.60 | 0.94 – 13.75 | 0.061 | 2.24 | 0.49 – 10.19 | 0.294 | 2.89 | 0.89 – 9.31 | 0.076 |
| Immature|Mature | 4.40 | 1.14 – 16.94 | 0.032 | 2.75 | 0.60 – 12.52 | 0.190 | 3.52 | 1.08 – 11.43 | 0.037 |
| Mature|2nd Cauliflower | 95.60 | 19.81 – 461.41 | <0.001 | 61.98 | 11.36 – 338.24 | <0.001 | 63.92 | 16.10 – 253.77 | <0.001 |
| Temperature [linear] | 1.53 | 0.87 – 2.73 | 0.145 | 9.54 | 0.99 – 104.48 | 0.057 | 1.47 | 0.85 – 2.59 | 0.174 |
| Temperature [quadratic] | 0.42 | 0.23 – 0.76 | 0.006 | 0.93 | 0.09 – 9.70 | 0.949 | 0.45 | 0.25 – 0.81 | 0.009 |
| Isolate [C1] | 1.58 | 0.59 – 4.28 | 0.366 | 1.49 | 0.55 – 4.03 | 0.432 | |||
| Isolate [C18] | 0.39 | 0.15 – 1.02 | 0.060 | 0.36 | 0.13 – 0.94 | 0.042 | |||
| Isolate [C24] | 0.70 | 0.26 – 1.83 | 0.464 | 0.62 | 0.23 – 1.66 | 0.346 | |||
| DeathAge | 1.11 | 1.07 – 1.14 | <0.001 | 1.10 | 1.06 – 1.14 | <0.001 | 1.09 | 1.06 – 1.12 | <0.001 |
|
DeathAge × Temperature [linear] |
0.96 | 0.92 – 1.01 | 0.107 | ||||||
|
DeathAge × Temperature [quadratic] |
0.98 | 0.93 – 1.03 | 0.517 | ||||||
| Observations | 162 | 162 | 162 | ||||||
| R2 Nagelkerke | 0.348 | 0.362 | 0.299 | ||||||
brant(modele14)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 16.1 18 0.59
## Temperature.L 4.04 3 0.26
## Temperature.Q 6.78 3 0.08
## IsolateC1 0.22 3 0.97
## IsolateC18 -461.22 3 1
## IsolateC24 2.63 3 0.45
## DeathAge 9.09 3 0.03
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
# GSimplified stages
FC1 = polr(stade_simplifie ~ Isolate + Temperature * DeathAge,
data = infected, Hess = TRUE)
model_selection = dredge(FC1)
model_selection
## Global model call: polr(formula = stade_simplifie ~ Isolate + Temperature * DeathAge,
## data = infected, Hess = TRUE)
## ---
## Model selection table
## (Int) DtA Isl Tmp DtA:Tmp df logLik AICc delta weight
## 8 + 0.09511 + + 8 -125.533 268.0 0.00 0.435
## 16 + 0.08670 + + + 10 -123.565 268.6 0.58 0.325
## 6 + 0.08108 + 5 -130.010 270.4 2.40 0.131
## 14 + 0.07302 + + 7 -128.273 271.3 3.27 0.085
## 4 + 0.08703 + 6 -131.054 274.7 6.64 0.016
## 2 + 0.07531 3 -134.837 275.8 7.82 0.009
## 5 + + 4 -148.833 305.9 37.91 0.000
## 1 + 2 -151.702 307.5 39.47 0.000
## 7 + + + 7 -148.297 311.3 43.31 0.000
## 3 + + 5 -151.140 312.7 44.66 0.000
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| stade_simplifie | stade_simplifie | |||||
|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p |
| Caul Grapes|Imm Mature | 3.04 | 0.79 – 11.66 | 0.105 | 1.72 | 0.38 – 7.84 | 0.483 |
| Imm Mature|2nd Caul | 71.72 | 15.03 – 342.18 | <0.001 | 43.91 | 8.21 – 234.87 | <0.001 |
| DeathAge | 1.10 | 1.07 – 1.14 | <0.001 | 1.09 | 1.06 – 1.13 | <0.001 |
| Isolate [C1] | 1.54 | 0.57 – 4.15 | 0.392 | 1.43 | 0.53 – 3.87 | 0.477 |
| Isolate [C18] | 0.40 | 0.15 – 1.03 | 0.063 | 0.36 | 0.13 – 0.94 | 0.042 |
| Isolate [C24] | 0.77 | 0.29 – 2.05 | 0.604 | 0.69 | 0.25 – 1.84 | 0.458 |
| Temperature [linear] | 1.57 | 0.89 – 2.82 | 0.126 | 14.72 | 1.41 – 164.42 | 0.028 |
| Temperature [quadratic] | 0.45 | 0.25 – 0.82 | 0.010 | 1.51 | 0.15 – 15.81 | 0.728 |
|
DeathAge × Temperature [linear] |
0.95 | 0.91 – 1.00 | 0.056 | |||
|
DeathAge × Temperature [quadratic] |
0.97 | 0.92 – 1.02 | 0.318 | |||
| Observations | 162 | 162 | ||||
| R2 Nagelkerke | 0.326 | 0.347 | ||||
best_mod = polr(stade_simplifie ~ Isolate + Temperature + DeathAge,
data = infected, Hess = TRUE)
tab_model(best_mod)
| stade_simplifie | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| Caul Grapes|Imm Mature | 3.04 | 0.79 – 11.66 | 0.105 |
| Imm Mature|2nd Caul | 71.72 | 15.03 – 342.18 | <0.001 |
| Isolate [C1] | 1.54 | 0.57 – 4.15 | 0.392 |
| Isolate [C18] | 0.40 | 0.15 – 1.03 | 0.063 |
| Isolate [C24] | 0.77 | 0.29 – 2.05 | 0.604 |
| Temperature [linear] | 1.57 | 0.89 – 2.82 | 0.126 |
| Temperature [quadratic] | 0.45 | 0.25 – 0.82 | 0.010 |
| DeathAge | 1.10 | 1.07 – 1.14 | <0.001 |
| Observations | 162 | ||
| R2 Nagelkerke | 0.326 | ||
brant(best_mod) # Not much better, but on death age and omnibus, same concluisons, so I kept the other one.
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 12.46 6 0.05
## IsolateC1 0.04 1 0.85
## IsolateC18 0.97 1 0.32
## IsolateC24 1.41 1 0.24
## Temperature.L 1.19 1 0.28
## Temperature.Q 0.72 1 0.4
## DeathAge 9.95 1 0
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
tab_model(modele14)
| stade_max | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| Cauliflower|Grapes | 3.24 | 0.85 – 12.35 | 0.084 |
| Grapes|Immature | 3.60 | 0.94 – 13.75 | 0.061 |
| Immature|Mature | 4.40 | 1.14 – 16.94 | 0.032 |
| Mature|2nd Cauliflower | 95.60 | 19.81 – 461.41 | <0.001 |
| Temperature [linear] | 1.53 | 0.87 – 2.73 | 0.145 |
| Temperature [quadratic] | 0.42 | 0.23 – 0.76 | 0.006 |
| Isolate [C1] | 1.58 | 0.59 – 4.28 | 0.366 |
| Isolate [C18] | 0.39 | 0.15 – 1.02 | 0.060 |
| Isolate [C24] | 0.70 | 0.26 – 1.83 | 0.464 |
| DeathAge | 1.11 | 1.07 – 1.14 | <0.001 |
| Observations | 162 | ||
| R2 Nagelkerke | 0.348 | ||
## `summarise()` has grouped output by 'Temperature', 'Isolate'. You can override
## using the `.groups` argument.
The quadratic effect of temperature (OR = 0.42, p = 0.006) indicates a
non-linear relationship, where parasite maturation is optimal at
intermediate temperatures (30°C) and decreases at both lower (20°C) and
higher (40°C) extremes. At extreme temperatures, the odds of reaching a
more advanced maturation stage decrease by 58%. There is also a non
significant trend of a less efficient maturation for C18, which is
consistent with the previous analysis for the second stage cauliflower
formation.
This is a supplementary analysis that we decided not to include as it would be outside the scope of our article.
FC1 = glmmTMB(data=data_mat, mature_load/100000 ~ Never_repro * Temperature + Isolate ,family= gaussian(), REML=T)
FC1R = glmmTMB(data=data_mat, mature_load/100000 ~ Never_repro * Temperature + Isolate + (1|Batch), family= gaussian(), REML=T)
AICc(FC1, FC1R)
## df AICc
## FC1 10 1424.867
## FC1R 11 NA
data_mat$Nb_offspring = !is.na(data_mat$Nb_offspring)
FC1 = glmmTMB(data=data_mat, mature_load/100000 ~ Never_repro * Temperature * Isolate * DeathAge,
family= gaussian(), na.action = na.fail)
model_selection = dredge(FC1)
head(model_selection)
## Global model call: glmmTMB(formula = mature_load/1e+05 ~ Never_repro * Temperature *
## Isolate * DeathAge, data = data_mat, family = gaussian(),
## na.action = na.fail, ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Nvr_rpr) cnd(DtA:Isl)
## 8 -48.740 + 1.9350 + 10.920
## 40 -35.850 + 1.6910 + -19.910
## 56 15.510 + 0.6365 + -28.080 +
## 4 -39.320 + 1.8480 +
## 24 -2.410 + 0.9804 + 9.321 +
## 20 9.883 + 0.8171 + +
## cnd(DtA:Nvr_rpr) df logLik AICc delta weight
## 8 7 -705.627 1426.1 0.00 0.303
## 40 0.6523 8 -704.766 1426.6 0.52 0.233
## 56 0.7908 11 -701.897 1427.8 1.71 0.129
## 4 6 -707.590 1427.8 1.72 0.128
## 24 10 -703.152 1427.9 1.87 0.119
## 20 9 -704.605 1428.5 2.47 0.088
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| mature_load/1e+05 | mature_load/1e+05 | mature_load/1e+05 | mature_load/1e+05 | mature_load/1e+05 | |||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| DeathAge | 1.93 | 1.45 – 2.42 | <0.001 | 1.69 | 1.08 – 2.30 | <0.001 | 0.64 | -0.48 – 1.75 | 0.264 | 1.85 | 1.36 – 2.34 | <0.001 | 0.98 | -0.06 – 2.02 | 0.065 |
| Isolate [C1] | 7.59 | -7.98 – 23.15 | 0.340 | 6.74 | -8.79 – 22.27 | 0.395 | -60.72 | -135.66 – 14.22 | 0.112 | 8.57 | -7.18 – 24.32 | 0.286 | -61.46 | -137.05 – 14.12 | 0.111 |
| Isolate [C18] | 28.90 | 14.05 – 43.74 | <0.001 | 28.18 | 13.38 – 42.98 | <0.001 | -46.60 | -112.61 – 19.41 | 0.166 | 27.62 | 12.62 – 42.61 | <0.001 | -37.89 | -103.59 – 27.81 | 0.258 |
| Isolate [C24] | 9.61 | -5.91 – 25.13 | 0.225 | 9.10 | -6.35 – 24.55 | 0.248 | -31.03 | -97.81 – 35.75 | 0.362 | 8.56 | -7.14 – 24.26 | 0.286 | -30.78 | -98.14 – 36.58 | 0.370 |
| Never repro | 10.92 | 0.19 – 21.65 | 0.046 | -19.91 | -67.05 – 27.23 | 0.408 | -28.08 | -75.34 – 19.18 | 0.244 | 9.32 | -1.34 – 19.98 | 0.087 | |||
| DeathAge × Never repro | 0.65 | -0.32 – 1.62 | 0.188 | 0.79 | -0.18 – 1.76 | 0.112 | |||||||||
| DeathAge × Isolate [C1] | 1.46 | -0.18 – 3.10 | 0.081 | 1.51 | -0.15 – 3.16 | 0.074 | |||||||||
| DeathAge × Isolate [C18] | 1.51 | 0.21 – 2.82 | 0.023 | 1.36 | 0.06 – 2.66 | 0.040 | |||||||||
| DeathAge × Isolate [C24] | 0.84 | -0.50 – 2.18 | 0.219 | 0.85 | -0.50 – 2.20 | 0.219 | |||||||||
| Observations | 145 | 145 | 145 | 145 | 145 | ||||||||||
| R2 / R2 adjusted | 0.395 / 0.368 | 0.402 / 0.371 | 0.425 / 0.382 | 0.378 / 0.356 | 0.415 / 0.376 | ||||||||||
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
## Family: gaussian ( identity )
## Formula: mature_load/1e+05 ~ DeathAge + Isolate + Never_repro + 1
## Data: data_mat
##
## AIC BIC logLik deviance df.resid
## 1425.3 1446.1 -705.6 1411.3 138
##
##
## Dispersion estimate for gaussian family (sigma^2): 987
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -48.739 13.903 -3.506 0.000456 ***
## DeathAge 1.935 0.249 7.770 7.87e-15 ***
## IsolateC1 7.586 7.944 0.955 0.339599
## IsolateC18 28.896 7.576 3.814 0.000137 ***
## IsolateC24 9.611 7.920 1.213 0.224975
## Never_repro 10.923 5.476 1.995 0.046064 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.348 0.352 0.484 0.78 0.292 0.672 0.328 0.38 0.696 0.468 0.104 0.048 0.164 0.24 0.468 0.428 0.596 0.376 0.268 0.412 ...
It seems that there is a possible slight positive impact of castration
on parasite load. However, the effective in the different categories are
really small.
## `summarise()` has grouped output by 'Never_repro', 'Isolate'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Never_repro', 'Isolate'. You can override
## using the `.groups` argument.
sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Ventura 13.2.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Paris
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] brant_0.3-0 ordinal_2023.12-4.1 car_3.1-3
## [4] carData_3.0-5 MASS_7.3-64 brms_2.22.0
## [7] Rcpp_1.0.14 coxme_2.2-22 bdsmatrix_1.3-7
## [10] glmmTMB_1.1.10 forcats_1.0.0 tidyr_1.3.1
## [13] dplyr_1.1.4 gridExtra_2.3 survminer_0.5.0
## [16] ggpubr_0.6.0 survival_3.8-3 MuMIn_1.48.4
## [19] sjPlot_2.8.17 DHARMa_0.4.7 patchwork_1.3.0
## [22] ggplot2_3.5.1 readr_2.1.5
##
## loaded via a namespace (and not attached):
## [1] tensorA_0.36.2.1 jsonlite_1.8.9 datawizard_1.0.2
## [4] magrittr_2.0.3 estimability_1.5.1 farver_2.1.2
## [7] nloptr_2.1.1 rmarkdown_2.29 vctrs_0.6.5
## [10] minqa_1.2.8 effectsize_1.0.0 rstatix_0.7.2
## [13] htmltools_0.5.8.1 distributional_0.5.0 broom_1.0.7
## [16] Formula_1.2-5 sjmisc_2.8.10 sass_0.4.9
## [19] bslib_0.9.0 plyr_1.8.9 emmeans_1.10.7
## [22] zoo_1.8-13 cachem_1.1.0 TMB_1.9.17
## [25] commonmark_1.9.5 mime_0.12 iterators_1.0.14
## [28] lifecycle_1.0.4 pkgconfig_2.0.3 gap_1.6
## [31] sjlabelled_1.2.0 Matrix_1.7-2 R6_2.5.1
## [34] fastmap_1.2.0 shiny_1.10.0 rbibutils_2.3
## [37] digest_0.6.37 numDeriv_2016.8-1.1 colorspace_2.1-1
## [40] qgam_1.3.4 labeling_0.4.3 km.ci_0.5-6
## [43] abind_1.4-8 mgcv_1.9-1 compiler_4.4.2
## [46] doParallel_1.0.17 bit64_4.6.0-1 withr_3.0.2
## [49] backports_1.5.0 performance_0.13.0 ggsignif_0.6.4
## [52] sjstats_0.19.0 ucminf_1.2.2 loo_2.8.0
## [55] tools_4.4.2 httpuv_1.6.15 glue_1.8.0
## [58] promises_1.3.2 nlme_3.1-167 gridtext_0.1.5
## [61] grid_4.4.2 checkmate_2.3.2 generics_0.1.3
## [64] gtable_0.3.6 tzdb_0.4.0 KMsurv_0.1-5
## [67] data.table_1.16.4 hms_1.1.3 xml2_1.3.6
## [70] foreach_1.5.2 pillar_1.10.1 markdown_2.0
## [73] stringr_1.5.1 vroom_1.6.5 later_1.4.1
## [76] posterior_1.6.1 splines_4.4.2 ggtext_0.1.2
## [79] lattice_0.22-6 bit_4.5.0.1 tidyselect_1.2.1
## [82] knitr_1.49 reformulas_0.4.0 litedown_0.6
## [85] stats4_4.4.2 xfun_0.51 bridgesampling_1.1-2
## [88] matrixStats_1.5.0 stringi_1.8.4 yaml_2.3.10
## [91] boot_1.3-31 codetools_0.2-20 evaluate_1.0.3
## [94] tibble_3.2.1 cli_3.6.3 RcppParallel_5.1.10
## [97] xtable_1.8-4 parameters_0.24.2 Rdpack_2.6.2
## [100] munsell_0.5.1 jquerylib_0.1.4 survMisc_0.5.6
## [103] ggeffects_2.2.1 coda_0.19-4.1 parallel_4.4.2
## [106] rstantools_2.4.0 bayestestR_0.15.2 gap.datasets_0.0.6
## [109] bayesplot_1.11.1 Brobdingnag_1.2-9 lme4_1.1-36
## [112] mvtnorm_1.3-3 scales_1.3.0 insight_1.1.0
## [115] purrr_1.0.4 crayon_1.5.3 rlang_1.1.5